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ABSTRACT 


Low  frequency  acoustic  propagation  in  shallow  water  is  examined  from  a  normal 
mode  context.  By  modelling  the  far  field  pressure  field  as  a  modal  sum,  propagating  mode 
characteristics  of  wavenumber,  initial  phase,  attenuation  and  amplitude  may  be  estimated 
using  a  high  resolution  parameter  modeling  technique.  The  advantages  of  such  an 
algorithm  are  the  resolution  of  closely  spaced  modes  in  a  range  independent  environment 
and  the  ability  to  analyze  range  dependent  waveguides. 

This  thesis  presents  the  application  of  a  Prony  algoiiilim  to  the  shallow  water 
environment.  The  algorithm  operates  directly  on  the  signal  matrix.  Synthetically 
generated,  range  independent  pressure  fields  are  used  to  analyze  the  technique's 
performance  and  to  observe  its  sensitivity  to  variations  in  model  specifications.  Noise  is 
added  to  determine  the  threshold  of  acceptable  pierformance.  As  a  consequence  of  field  data 
tests,  further  enhancements  to  the  algorithm  are  suggested. 

Range  dependent  performance  is  evaluated  on  a  coastal  wedge  example  and 
geoacoustic  parameter  shift  example. 

Thesis  advisor:  George  V.  Frisk,  Associate  Scientist,  Woods  Hole  Oceanographic 
Institution 
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Chapter  1 
Introduction 


1.1  Introduction 

In  this  thesis,  we  shall  investigate  the  application  of  a  signal  processing  technique  to 
acoustic  propagation  in  a  shallow  water  waveguide.  As  is  the  case  in  many  disciplines, 
advances  in  the  ocean  acoustics  community  tend  to  occur  in  incremental  steps.  Often,  the 
approaches  used  in  acoustic  problem  formulation  or  analysis  are  successfully  applied  in  other 
research  efforts.  Examples  include  ray  theory  (from  optics),  method  of  images 
(electromagnetic  theory )  and  fast  field  programs  (digital  signal  processing).  We  shall  examine 
the  application  of  Prony's  method  to  the  problem  of  resolving  normal  modes  in  CW  data 
obtained  on  a  horizontal  array  in  shallow  water.  We  expect  to  formulate  a  technique  which 
requires  the  researcher  to  incorporate  knowledge  of  the  acoustic  propagation  into  a  signal 
modelling  problem.  In  return  for  this  bounding  or  constraint  of  the  problem,  we  anticipate  a 
gain  in  two  areas.  First,  in  the  range  independent  waveguide,  Prony's  method  should  allow 
the  resolution  of  closely  spaced  propagating  modes.  Second,  in  a  range  dependent  waveguide, 
the  short  apertures  used  by  the  high  resolution  techniques  may  permit  estimation  of  waveguide 
parameters  by  an  adiabatic  assumption. 

To  obtain  such  high  resolution  using  conventional  discrete  Hankel  transforms  and  FFT 
beamforming  techniques  requires  a  large  aperture  since  the  resolution  is  inversely  proportional 
to  the  length  of  the  array[l].  In  theoretical  or  computer  generated  fields,  this  lengtli  of  data 
may  be  easily  obtained.  Shallow  water  data  rarely  meets  this  criterion  since  ocean  waveguides 
which  are  invariant  over  an  interval  of  kilometers  (which  is  needed  to  obtain  the  desired 
resolution)  are  not  generally  found  in  practice. 

Applications  of  Prony's  method  may  be  found  in  a  variety  of  disciplines.  Research 
examples  include  seismic  exploration|2),  acoustic  echo  reductionf3],  subsurface  radar[41, 

-.S- 
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beamforming[5],  and  structural  analysis[61 .  The  specific  use  of  Prony's  method  in  a  shallow 
water  waveguide  modal  context  has  been  explored  in  parallel  with  this  thesis  development  by 
Shang,  et  al.[7].  The  Shang  study  considers  the  algorithm  in  a  scheme  to  localize  a  source  and 
provide  some  waveguide  characterization.  Their  characterization  approach  differs  from  this 
thesis  in  order  selection,  filter  coefficient  determination  and  attenuation  parameter  estimation. 
Such  differences  in  implementation  emphasizes  the  flexible  nature  of  algorithm  structure  within 
the  framework  of  Prony's  method. 

The  motivation  for  exploring  this  approach  is  to  assist  in  the  effort  to  solve  the  inverse 
problem  of  determining  geoacoustic  parameters.  Specifically,  the  geoacoustic  parameter 
inverse  technique  requires  the  set  of  horizontal  wavenumbers  of  the  propagating  modes  as 
input  data[8-l  1].  The  signal  processing  method  used  provides  a  means  to  identify  the 
horizontal  wavenumber  of  an  acoustic  field.  Ultimately,  we  envision  measurements  of  the 
pressure  field  in  the  water  column  yielding  a  set  of  wavenumbers  which  may  then  be  used  to 
infer  bottom  properties. 

There  are  two  sets  of  field  data  which  have  been  collected  to  date.  These  provide  a  test 
envuronment  of  the  algorithm  in  the  real  world.  The  first  set  was  collected  in  May  of  1984  off 
the  coast  of  Nantucket  Island,  MA  at  140  Hz  and  220  Hz[10]  and  the  other  set  was  obtained  in 
September  of  1985  off  the  coast  of  Corpus  Christi,  TX  at  50  Hz  and  140  Hz[12].  These  sites 
were  selected  since  the  bottoms  are  reasonably  flat  and  other  studies  of  bottom  properties  are 
available  for  comparison  of  the  geoacoustic  parameters.  The  experimental  setup  uses  a  small 
vessel  which  tows  a  two  frequency  CW  source  at  a  fixed  depth[8].  The  source  is  towed  away 
from  two  fixed  moored  receivers  and  the  horizontal  source-receiver  range  is  monitored  via  a 
radar  tracking  system.  The  data  collection  technique  consists  of  drifting  away  from  the 
receivers.  The  drift  rate  must  be  low  (<  0.5  knots)  to  provide  adequate  sampling  (to  avoid 
aliasing)  of  the  pressure  field.  The  hydrophones  are  part  of  a  larger  BODIS  assembly  which 
removes  harmonic  time  dependence  of  the  pressure  data  by  quadrature  demodulation. 
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The  effect  of  the  drift  scheme  and  quadrature  demodulation  is  to  model  the  data  set  as  being 
collected  on  a  synthetic  array.  The  physical  setup  is  summarized  in  figure  1.1.1. 


Fig  1.1.1  Experimental  setup  for  data  collection 


7 


1.2  Outline 


Chapter  2  provides  a  basic  review  of  normal  mode  theory  in  an  acoustic  waveguide  by 
developing  the  mathematical  description  in  terms  of  a  Stumi  Liouville  problem.  The  second 
section  provides  a  derivation  of  the  Prony  model  and  explores  the  nuances  of  the  algorithm 
used.  The  last  part  of  the  chapter  casts  the  normal  mode  environment  in  terms  of  the  model  and 
illustrates  its  validity.  Two  examples  are  provided:  the  first  involves  the  application  of  the 
method  to  synthetic  pressure  field  data,  while  the  second  demonstrates  the  ability  of  the 
technique  to  extract  mode  shapes  through  use  of  a  vertical  array. 

In  Chapter  3,  the  depth  dependent  Green's  function  is  reviewed  and  the  Prony  energy 
spectral  density  (ESD)  is  introduced.  The  ESD  is  used  as  a  tool  to  transform  all  parameter 
estimates  into  a  simple  graphical  display.  As  such,  it  acts  as  a  tool  to  aid  in  waveguide 
analysis.  Additional  tools  are  developed  and  examined,  and  the  algorithm  is  tested  using  these 
tools  on  two  synthetically  generated  sets  of  data.  The  effect  of  changes  in  input  variables  on 
parameter  estimates  is  examined  and  ranked  according  to  sensitivity.  This  ranking  leads  to  the 
development  of  a  set  of  empirically  derived  guidelines  to  specify  a  model  order,  ap)erture  size 
and  averaging.  The  chapter  concludes  with  a  performance  evaluation  on  field  data. 

Range  dependent  performance  is  addressed  in  Chapter  4.  Two  examples  of  a  range- 
dependent  environment  are  provided.  One  waveguide  consists  of  an  upslope  propagation  in  a 
coastal  wedge  scenario.  The  second  contains  a  step  change  in  bottom  geoacoustic  properties. 

Chapter  5  provides  a  "wish  list"  of  enhancements  and  considerations  in  further  versions 
of  this  high  resolution  scheme.  Topping  the  list  is  the  issue  of  more  robust  performance  in 
noisy  environments.  An  example  of  a  bandpa.ss  filter  scheme  to  address  this  issue  is  also 
provided. 


Chapter  2 
Basic  Theory 


2.1  Normal  Mode  Theory  Review 

A  review  of  basic  normal  mode  theory  is  important  to  emphasize  assumptions  made 
in  the  development  of  the  description  ol  the  acoustic  field  and  in  the  application  of  the 
Prony  model  to  this  environment  An  advantage  of  the  normal  mode  approach  is  that  it 
allows  us  to  build  on  previous  results  as  the  boundary  conditions  become  more  complex. 

It  will  be  shown  that  this  method  makes  the  problem  tractable  by  casting  the  linear,  second 
order  differential  equation  arising  from  the  boundary  value  problem  in  terms  of  a  Sturm- 
Liouville  problem.  This  allows  us  to  apply  the  rich  existing  mathematical  theory  for  the 
subject. 

In  this  thesis,  sound  propagation  in  shallow  water  will  be  treated  as  a  field 
propagating  within  a  waveguide  constrained  by  the  surface  and  the  bottom.  The  field 
distribution  within  the  waveguide  is  affected  by  the  boundary  characteristics  so  that  the 
local  modes  act  as  a  sampling  mechanism  for  the  properties  of  the  boundaries.  It  is 
assumed  that  the  top  boundary  condition  remains  constant  (in  this  case,  a  pressure  release 
surface)  while  the  mode  wavenumbers  are  affected  by  changes  in  the  dimension  of  the 
waveguide  (bathymetry)  and  bottom  properties  (such  as  sound  speed,  density  and 
attenuation)  (sec  fig  2.1.1).  The  water  column  is  the  region  of  interest  for  measurement 
purposes  since  the  pressure  field  may  be  easily  obtained  by  real  towed  or  synthetic  aperture 
arrays.  In  this  section,  we  will  assume  that  the  wavegu>.  e  is  locally  range-independent. 
After  presenting  some  general  normal  mode  theory,  this  section  will  use  hard  bottom  and 
Pekeris  waveguidc.s  to  emphasi/.c  salient  points 
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Fig  2.1.1.  Shallow  Water  Waveguide  Model 


In  the  ocean  acoustics  waveguide,  the  governing  equation  for  ilie  acoustic  pressure, 


P,  is  the  time  dependent  wave  cquation[13]: 
1  82 
c'-tD  3t2 


(2.1.1) 


where  F(r,t)  is  a  source  function.  For  the  case  of  interest,  we  will  consider  a  harmonic 
point  source  function.  We  therefore  let  PCr.t)  =  p(r)e'J*^^  and  F(r,t)=f(r)c  which 
transforms  the  wave  equation  to  the  inhomogeneous  Helmholtz  equation: 

(V2  +  k2(r))p(i)  =  -47C  f(r)  .  k(D  =  ^  (2. 1 .2) 

The  point  source  is  modeled  as  an  impulse  function  of  strenglli  S,  ie: 

(V2  +  k2(r))p(r)  = -4a:  S  6(r  -  ro)  (2.1.3) 


In  cylindrical  coordinates,  assuming  horizontal  stratification  k(r)  =  k(z)  and  a  with  (})(r)  = 


p(x)  and  a  source,  located  at  (ro  7.<)_B()),  the  Heimliolt/  equation  is]  14]: 
1  8  8p  I  82p  f)?.p 

- (r — )  4  — — -  +  — -+  k-(/.)  p  -471 - o(/ 

^8r  8r  r2ao?.  o/.2  '  ^ 


471  7,-d  sai-Op)  (2.;. 4) 
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Assuming  cylindrical  symmetry,  p(r)  =  p(r,z),  and  ro=0,the  Helmholtz  equation  may  be 

integrated  with  respect  to  0  to  remove  angular  dependency: 

2k 

J'l-p-  ^(r  pi  ^5(z-zo)  (2.1 

0 


which  leads  to: 


1  9  ^  3p.  02p  -  5(r)  _ 

p  =  -2— 


(2.1.6) 


The  pressure  field  is  constrained  by  the  following  boundary  conditions: 

•  p(r,0)  =  0  (pressure  release  surface) 

•  a  bottom  impedence  boundary  condition 


and  a  Sommerfeld  radiation  condition  (which  specifies  energy  from  the  source  as 
propagating  outward)!  14]. 

The  two  dimensional  Green's  function  can  be  expanded  in  a  complete  orthonormal 
set  of  the  eigenfunctions  of  depth,  z[15,16].  The  method  of  solving  an  inhomogeneous 
Sturm  Liouville  equation  (see  Appendix  A)  may  be  used  as  follows.  First,  the 
eigenfunctions,  Un(z),  are  found  by  solving  the  homogeneous  Helmholtz  equation  by 


separation  of  variables.  This  yields  an  equation  for  the  Un(z): 
d^u  2  2 

^j-Y+  fk^  -  k^ju  =  0  and  k^  is  the  separation  constant 


(2.1.7) 


The  solution  for  p  is  then  assumed  to  have  the  form  p(r,z)  =  ZRn(r)un(z)  and  the  assumed 

n=l 

solution  is  substituted  into  equation  (2. 1 .6),  the  inhomogeneous  Helmholtz  equation. 

Then,  the  orthogonality  and  completeness  characteristics  of  the  eigenfunctions,  Un(z),  ;ue 
used  to  find  the  coefficients  Rn(r).  For  an  arbitrary  k(z),  the  solution  may  lx:  expressed  as, 

p(r,7.)  Jr.  ^^u^^(/.())  u„(7.)  Il*''l(k„r)  +  l(r)  (2. 1 .8) 

n  -1 
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in  which  H^Q^(knr)  is  a  Hankel  function  of  the  first  type  and  is  generated  by  the  solution  to 


equation  (2.1.7).  While  the  discrete  sum  corresponds  to  the  trapped  modes,  the  continuum 
contribution,  I(r),  consists  of  branch  line  integrals  and  "improper"  modes  (depending  on 
the  branch  cut  selection)[17-20]. 

Restricting  our  attention  to  the  hard  bottom  case  requires  imposing  a  bottom 
3p  I 

boundary  condition  of —  Iz  =  h  =  0.  The  hard  bottom  isovelocity  case  yields  an  analytic 

9z 

solution  for  Un(z)[14]: 

oo 

p(r,z)  =  jtc  ^  sin(7nZ0)  sin(7nz)  q  ^ 

n=l 


=  vertical  wavenumber 
(n  -  It 
h 

Although  a  perfectly  reflecting  (hard)  bottom  is  not  found  in  nature,  this  simple 
model  allows  identification  of  three  important  normal  mode  characteristics  which  may  be 
used  in  more  complex  models.  First,  as  the  range  from  the  source  increases,  the  Hankel 
function  may  be  replaced  by  an  asymptotic  approximation[21 ): 

H%knr)  =-v/-7^ei(knr-J)  (2.1.10) 

\  ttknr 

Second,  a  propagating  mode  may  be  viewed  as  a  cylindrically  spreading,  outgoing  wave 
with  a  vertical  shape  determined  by  the  mode  eigenfunctions,  Un(z)[221.  'Fliird,  the 
placement  of  source  and  receiver  directly  affects  mode  amplitudes.  In  equations  (2. 1 .8) 
and  (2.1.9),  we  can  see  the  product  of  the  source  and  receiver  eigenfunctions  determine  the 
mode  excitation  available  at  the  receiver. 

A  more  complex  model  introduced  by  Pekeris  in  1948  consists  of  approximating 
the  bottom  by  a  half  space  with  a  constant  density  and  .sound  speodl  23),  The  resulting 
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eigenvalue  spectrum  consists  of  two  seperate  regions;  one  region  of  the  spectrum  has  a 
discrete  spectrum  and  the  other  region  has  a  continuous  spectrum.  This  change  in  spectrum 
characteristics  is  due  to  the  bottom  boundary  defined  as  an  impedence  condition; 

po3z  Pb  3z 
P  L  =  h  =  Pb  lz  =  h 

The  horizontal  wavenumber  which  separates  the  discrete  and  continuous  regions  is  known 
as  the  cutoff  wavenumber  is  delineated  by  the  wavenumber  of  the  bottom  half  space,  which 
has  sound  speed  Cb  and  density  pb-  This  discrete/continuous  split  spectrum  is  found  in 
more  complex  models  in  which  the  horizontally  stratified  layers  that  make  up  the  bottom 
profile  are  usually  terminated  at  an  arbitrary  depth  by  an  isovelocity  halfspace.  The  cutoff 
wavenumber  is  set  by  this  layer's  wavenumber [24]. 

An  intutive  approach  to  cutoff  is  formulated  by  alternately  posing  the  modal 
description  as  a  superposition  of  up  and  down  going  plane  waves(since  Un(z)  =  sin(ynz)  = 

gjTfnz  .  g-jTnz 

- 2j - •)[  13,25].  This  expression  is  well  suited  to  identifying  modes  which  will 

propagate. 


Fig  2.1.2  -  Plane  Waves  Incident  on  Bottom 


As  shown  in  figure  2.1.2,  the  downgoing  plane  wave  is  incident  on  the  bottom  at  an  angle 
6n,  which  is  defined  by  0n  =  tan'^f— )  .  As  the  mode  number,  n,  increases,  the  angle  0n 

lYn  J 

decreases  and  so  the  inclination  of  the  plane  waves  becomes  more  nearly  vertical.  As  kn 
approaches  zero  and  becomes  imaginary,  the  mode  changes  from  a  propagating  mode  to  an 
exponentially  decaying,  inhomogeneous  wave.  This  situation  in  which  the  sound  energy  is 
present  as  a  heavily  attenuated  field  is  known  as  cutoff;  the  angle  at  which  this  occurs  is 
known  as  the  critical  angle,  0c[26].  An  incident  plane  wave  associated  with  a  propagating 
mode  will  experience  total  reflection  at  the  bottom;  in  the  sediment,  the  sound  pressure  is 
exponentially  damped  (Fig  2. 1 .3).This  "impedance  condition"  dictates  the  existence  of  a 
mechanism  to  account  for  energy  "leaking"  into  the  bottom  since  the  plane  waves  (the 
modes  are  being  modelled  as  a  superposition  of  up  and  downgoing  plane  waves)  incident 
on  the  bottom  are  no  longer  being  perfectly  reflected. 

By  decomposing  the  depth  eigenfunctions  into  plane  waves ,  the  plane  wave 
incident  angles  which  lead  to  propagating  modes  may  be  divided  into  two  regions[27]: 

•  0  >  0^  :  Region  of  perfect  reflection  resulting  in  a  discrete  set  of  trapped  modes  in 
water  column  and  exponential  decay  in  the  bottom. 

•  0  <  0r  :  Region  known  as  the  continuum  which  is  propagation  region  where  leaky 
(or  virtual)  modes  can  exist.  These  heavily  attenuated  modes  allow  energy  to  leak 
into  the  bottom;  some  of  this  energy  may  be  directed  back  into  the  waveguide  by 
the  bottom's  velocity/density  profile. 
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will  be  high  attenuation  with  respect  to  range.  While  the  water  column  mode  shapes 
won't  change  when  attenuation  is  taken  into  account,  different  modes  will  have  different 
attenuation  rates.  This  is  evident  if  the  plane  wave  decomposition  of  the  modes  approach 
is  used.  The  higher  order  modes  are  incident  on  the  bottom  with  steeper  angles  (closer  to 
normal  incidence).  These  modes  undergo  more  reflections  for  a  given  horizontal  range 
than  modes  which  are  incident  closer  to  grazing. 

The  cursory  review  of  normal  mode  theory  of  this  section  was  meant  to  empha.size 
key  aspects  of  the  shallow  water  waveguide  problem.  The  hard  bottom  waveguide 
example  demonstrated  the  modal  sum  form.  The  asymptotic  behavior  of  the  Hankel 
function  was  stated  as  well  as  characteristics  of  the  propagating  energy [27].  The 
substitution  of  an  impedance  condition  for  the  formerly  hard  bottom  results  in  two  distinct 
regions  in  the  spectrum.  In  one,  the  spectrum  is  discrete  and  in  the  other,  it  is  continuous. 
The  distinction  between  the  discrete  and  continuous  spectrum  may  be  set  by  the 
wavenumber  of  a  "basement"  isovelocity  halfspace.  The  general  pressure  field  description 
consists  of  a  modal  sum  and  continuum  as  in: 

oo 

p(r,z)  =  jti  ^u*(zo)  Un(z)  (knr)  +  I(r)  (2.1.1 2) 

n=l 


Attenuation  effects  are  incorporated  by  a  small  imaginary  term  in  the  horizontal 
wavenumber,  kn- 

If  the  asymptotic  form  of  the  Hankel  function  is  used  (equation  2.1.10)  and  the  attenuation 
term  is  included,  the  pressure  field  description  becomes: 


p(r,z)  -jtc 


P 

^^^aqU*(zo)  Uq(z) 

q=l 


V 


—  ei  -  7)  +  aq  r 
Ttr  ^ 


(2.1.13) 
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2.2  Prony’s  Method 

Irony's  method  is  a  parameter  estimation  technique  in  which  the  model  parameters 
are  varied  to  fit  the  observed  data.  Parameter  estimation  approaches  involve  the  use  of  a 
priori  knowledge  in  the  intelligent  selection  of  an  appropriate  model[29].  In  Prony's 
method,  the  signal  is  considered  to  be  composed  of  a  linear  combination  of  damped 
complex  exponentials[30].  Prony's  method  is  by  no  means  the  only  model  which  may  be 
applied  to  the  ocean  waveguide  (for  example  Pisarenko  or  autoregressive  moving  average 
(ARMA)  modelling  might  also  be  used)[31].  The  use  of  Prony's  algorithm  was  driven  by 
two  factors.  First,  the  method  has  the  advantage  of  requiring  short  data  lengths  ( a  small 
range  aperture)  to  yield  high  resolution  wavenumber  estimates.  Second,  the  modal 
structure  of  far  field  propagation  in  the  shallow  water  waveguide  can  fit  the  Prony  model 
very  well.  These  two  characteristics  lend  credibility  to  the  application  of  the  technique  to 
shallow  water  waveguide  propagation. 

In  1795,  Gaspard  Riche,  Baron  de  Prony,  proposed  an  interpolation  scheme  in 
which  a  deterministic  model  was  assumed  and  the  equally  spaced  data  was  used  to  fit  this 
deterministic  model[31].  The  method  consisted  of  an  exact  fit  of  the  data  points  to 
exponentials;  the  evolution  of  the  algorithm  since  then  has  been  significant.  The  insight  of 
the  solution  method  and  the  ensuing  three  step  process  has  endured  although  the  expanded 
algorithms  additionally  address  issues  such  as  stability,  robusmess  in  noise  and  least 
square  fits. 

The  model  used  by  the  algorithm  is  a  weighted  sum  of  complex  damped 
exponentials.  Consideration  of  the  exactly  determined  case  allows  identification  of  the 
steps  used  in  the  algorithm;  the  development  of  the  extended  Prony  method(least  squares 
fit  for  an  overdetermined  system)  is  an  enhancement  of  this  basic  procedure. 
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(2.2.1) 


The  observed  data,  y[n],  is  assumed  to  fit  the  model: 


with 


p  =  model  order 
T  =  sampling  range 
Aq  =  amplitude 
aq=  damping  factor 
kq  =  wavenumber 
0q  =  initial  phase 
Regrouping  terms: 


q=Q 


where 

Vq  =  Aqexp(j0q) 

Zq  =  exp[(aq+jkq)T] 

For  p  data  points,  the  system  may  be  expressed  in  matrix  form: 


(2.2.2) 


or 


/  y[0l  \ 
yH] 


JWJ 


Vy[p-ll/ 


(2.2.3) 


7.y_=X 

In  order  to  solve  for  the  complex  quantities  Vq  and  Zq,  we  will  decouple  equation  (  2.2.3) 
by  solving  for  the  Zq's.  Tlie  insight  offered  by  Prony  in  this  scheme  was  that  one  approach 
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for  finding  the  Zq's  is  the  solution  of  a  homogeneous  linear  constant  coefficient  difference 
equation  (LCCDE). 

The  first  two  steps  of  Prony's  method  solve  for  the  zq's  as  follows; 

Consider  the  observed  signal,y[n],  to  be  the  output  of  an  all  pole  filter  driven  by  an  unit 
impulse  function  (Fig  2.2.1): 


Fig  2.2.1  Filter  Model  of  Observed  Signal 


For  this  model: 


H(z) 


1  +  ^aqZ-q  f][(l-ZrZ-l) 
0=1  r=fl 


0=1 
,Y(Z) 


Since  H(z)  is  defined  as  cross  multiplying  equation  (2.2.4)  yields; 

Y(z)  1  +  i^aqz-q  =X(z) 

-  q=l  - 

Defining  ao  =  1  and  taking  the  inverse  Z  transform  results  in  an  LCCDE; 

P 

x[n]  =  Z  aq  y[n-q] 

0=0 

The  homogeneous  portion  of  the  LCCDE  is: 

i  aq  y[n-q]  =  0 
q=0 

Expanding  by  one  term  (recall  ao  =  1 ): 

P 

y[n]  +  Z  aqy[n-q]  =0 
0=1 


(2.2.4) 


(2.2.5) 


(2.2.6) 


(2.2.7) 


(2.2.8) 
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This  yields  a  set  of  equations  over  the  observed  data  range  0  <  n  <  2p-l: 


y[p-i] 

y[p-2] 

-  y[0] 

/ai\ 

/  y[pl  \ 

y[p] 

y[p-il 

...  y[i] 

a2 

=  - 

y[p+i] 

y[2p-2] 

y[2p-3] 

...  y[p-i] 

J 

Vapy 

^[2p-l]/ 

(2.2.9) 


The  first  step  of  Prony's  method  involves  obtaining  the  coefficients  aq's  by  solving 
equation  (2.2.9)  using  the  2p  complex  data  points  y[n].  Notice  that  the  aq's  were  defined 
in  equation  (2.2.4)  as  the  coefficients  of  a  polynomial  which  had  the  Zq's  as  the  roots. 
Obtaining  the  roots  of  the  equationis  the  second  step  of  the  process.  Now  that  the  Zq's  are 
available,  the  third  and  final  step  of  the  process  is  to  solve  the  exactly  determined  sytem 
defined  by  equation  (2.2.3).  The  nonlinear  aspect  of  the  problem  has  been  isolated  in  the 
second  step.  An  attempt  to  solve  equation  (2.2.1)  directly  by  an  error  minimization  scheme 
results  in  a  nonlinear  set  of  equations  which  must  be  solved  by  Newton's  method  or  some 
other  iterative  approach. 

In  the  last  section,  we  identified  the  three  steps  used  to  fit  observed  data  to  a  model 
defined  by  an  exactly  determined  set  of  equations.  In  practice,  the  amount  of  data  typically 
exceeds  the  model  order,  which  allows  extension  of  the  procedure  to  encompass  a  least 
squares  estimation.  The  advantage  of  the  least  squares  technique  is  that  the  issues  of  noise 
and  stability  may  be  approached  using  the  rich  theoretical  material  available  in  linear 
prediction,  lattice,  and  autoregressive  (AR)  filter  design[32].  This  section  will  take  a  closer 
look  at  the  three  steps  in  Prony's  method  and  identify  the  algorithm  used  in  our  work. 

The  first  step  contains  the  greatest  variety  of  approaches.  The  identification  of  tlie 
"best"  polynomial  coefficients  is  complicated  by  the  overdetermined  situation;  instead  of 
solving  equation  (2.2.9)  directly,  an  error  criterion  must  be  minimized.  There  are  three 
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general  ways  in  which  the  error  criterion  is  minimized,  one  "direct"  and  two  "indirect". 

The  use  of  "indirect"  and  "direct"  qualifiers  is  an  identification  of  how  the  coefficient  vector 
is  determined.  In  the  "direct"  method,  operations  take  place  directly  on  the  signal  matrix 
while  the  "indirect"  procedures  use  the  observed  data  to  generate  an  exactly  determined  set 
of  linear  equations.  The  "indirect"  methods  are  the  techniques  most  frequently  used  and 
will  be  discussed  first. 

The  "indirect"  methods  will  be  presented  here  from  the  viewpoint  of  linear 
prediction  theory[29,32-34].  A  simple  way  of  expressing  the  linear  prediction  philosophy 
is  that  the  p+1  output,  ie  y[p+l],  may  be  predicted  by  using  a  linear  combination  of  the  last 
p  outputs,  ie.  y[p],  y[p-l],...y[l].  By  using  our  all  pole  model  in  Figure  2.2.1,  we  can 
make  use  of  its  LCCDE  (equation  (2.2.6))  and  replace  the  input,  x[n],  by  5[n]. 


Since  ao  =  1,  LCCDE  may  be  rewritten  as: 

y[nl  =  -  I  aq  y[n-q]  +  6[n]  (2.2. 10) 

q=l 

Denoting  as  the  predicted  estimate  of  y[p+l]: 

t3  =  -  Zaqy[n-q]  (2.2.11) 

The  error,e(n],  is  now  defined  as  the  difference  between  the  observed  and  predicted  value: 

e[n]  =  y[n]-fr  (2.2.12) 

The  coefficients  aq  which  minimize  the  energy  in  the  error  signal  are  found  by  the  "normal 
equations".  In  matrix  form,  the  quantity  to  be  minimized  is: 

e  =  £'r£  (2.2.13) 

The  error  equations  of  equation  (2.2.12)  may  be  gathered  in  matrix  form: 

y-Ca  =  e  (2.2.14) 

Expanding  (2.2.13)  and  minimizing  yields  (he  following  relationship: 

C"Ca  =  CHy  (2.2.15) 
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which  are  known  as  the  "normal  equations".  These  equations  exploit  the  orthogonality  of 
the  error  vector  to  the  basis  set  formed  by  the  columns  of  the  signal  matrix.  The  range  of 
values  over  which  the  error  is  minimized  differentiates  the  two  "indirect"  methods,  which 
are  known  as  the  autocorrelation  and  autocovariance  methods[33]. 

The  autocorrelation  method  minimizes  the  error  over  an  infinite  duration  signal. 

The  method  assumes  a  stationary  process  and  yields  a  matrix  system  which  has  a  Toeplitz 
structure.  This  symmetric,  positive  semidefinite  characteristic  allows  the  use  of  Levinson's 
recursions  for  an  efficient  solution  algorithm.  The  main  advantage  of  the  autocorrelation 
method  is  that  the  ensuing  filter  is  theoretically  guaranteed  to  be  stable  (all  of  the  poles 
within  the  unit  circle  of  the  z  plane).  Care  must  be  taken  in  the  numerical  implementation  of 
the  technique  to  avoid  accumulated  roundoff  errors  from  making  the  autocorrelation  matrix 
ill  conditioned. 

The  disadvantages  stem  from  the  minimization  of  the  error  over  an  infinite  interval. 
In  an  all  pole  model,  the  impulse  response  will  be  infinite  in  duration.  The  finite  amount  of 
data  available  to  the  user  is  an  implicit  windowing  of  an  infinite  duration  signal.  This 
windowing  effect  changes  the  autocorrelation  coefficients  in  the  matrix,  forcing  an  estimate 
(rather  than  determination)  of  the  autocorrelation  coefficients.  The  result  of  spreading  the 
error  over  an  infinite  interval  is  that  the  model  generated  from  all  pole  data  will  not  match 
the  actual  system.  Zero  padding  and  application  of  windows  to  the  available  data  minimize 
the  effects  of  finite  data  length  but  these  techniques  may  be  hazardous  in  a  high  resolution 
spectrum  analyzer  situation. 

The  more  common  "indirect"  approach  in  use  today  is  the  autocovariance 
method(3,35,36].  The  error  is  minimized  over  the  finite  length  of  data.  In  stochastic 
theory  this  equates  to  the  nonstationary  case  modelled  as  locally  stationary.  The  matrix 
system  which  is  solved  for  the  aq’s  is  positive  semi-definite  but  not  Toeplit/..  Althougli 
Levinson's  recursions  cannot  be  used  to  solve  the  system,  tlicre  ;ire  algorithms,  most 
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notably  Marple's,  which  address  efficient  solution  techniques  for  this  method.  By 
minimizing  the  error  over  a  finite  interval,  the  autocovariance  method  will  match  all  pole 
data  to  the  generating  system.  As  the  data  length  increases,  the  covariance  method  will 
approach  the  autocorrelation  technique. 

A  drawback  to  the  covariance  method  is  the  lack  of  a  stability  theorem  to  guarantee 
the  filter  coefficients  will  describe  a  stable  filter.  This  drawback  was  addressed  by  Burg, 
who  developed  an  algorithm  which  constrained  the  problem  to  yield  a  stable  filter[37,38]. 
Specifically,  this  is  done  by  changing  the  error  minimization  problem.  In  addition  to  the 
error  defined  above,  called  the  forward  error,  a  new  error,  called  the  backwards  error,  is 
specified.  Burg's  rationale  was  that  the  stationary  signal  should  "look"  the  same  going 
forward  and  backwards  through  the  data  set  The  error  criteria  to  be  minimized  in  a  Burg 
algorithm  is  the  sum  of  the  forward  and  backward  errors.  Although  summing  the  error 
over  twice  as  many  points  is  advantageous  in  a  short  data  set,  blind  application  of  Burg’s 
algorithm  (also  known  as  the  modified  covariance  method)  is  dangerous.  Specifically,  the 
problem  lies  in  the  assumption  of  the  signal  appearing  the  same  regardless  of  the  direction 
of  data  set  traversal.  A  sinusoid  does  indeed  have  this  characteristic  but  the  presence  of 
damping  (placing  an  exponential  decaying  envelope  on  the  sinusoid)  requires  careful 
examination  of  the  physical  system  before  applying  the  Burg  forward -backward  error 
sum[2].  In  the  tradeoff  between  damping  and  stability,  the  decision  is  usually  made  to 
assume  an  undamped  system  in  return  for  the  assured  stability.  In  the  shallow  water 
waveguide,  constraining  the  pressure  field  to  consist  of  undamped  sinusoids  is  rea.sonable 
for  propagating  modes  but  not  for  the  leaky  or  virtual  modes  which  may  be  represented  as 
damped  sinusoids.  In  addition,  although  the  assumption  that  the  wavenumbers  are  real  (no 
damping)  is  acceptable  since  the  propagating  modes  have  small  damping  factors,  the  user 
must  tlicn  have  other  methods  available  to  estimate  nuxle  attenuation. 
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The  approach  explored  in  this  thesis  is  the  "direct"  method  in  which  the  least 
squares  operation  is  perforaied  directly  on  the  signal  matrix.  The  technique  is  hampered  by 
the  lack  of  a  stability  theorem  but  this  was  not  found  to  be  a  problem  in  practice.  In  tlie 
"direct"  approach,  a  QR  decomposition  is  used.  The  QR  decomposition  u  a  very  stable, 
fast  technique  which  uses  Householder  (orthogonal)  matrices  to  orthogonalize  the  original 
matrix[39,40].  The  use  of  Householder  matrices  results  in  a  decomposition  of  the  A 
matrix  into 

X  q  ~  Qp  X  p  X  (j  (2.2.16) 

in  which  tlie  first  q  columns  of  Q  form  an  orthonormal  basis  for  the  column  space  of  A  and 
the  last  p-q  columns  form  a  basis  for  the  left  nullspace  of  A  (ie,  the  last  p-q  columns  are 
perpendicular  to  the  first  q  columns).  The  first  n  rows  of  R  fonn  an  upper  triangular 
matrix;  the  columns  of  R  are  formed  by  successive  Householder  matrices  operating  on 
corresponding  columns  of  A. 

The  least  squares  problem  may  be  solved  by  considering; 

Ax  =  b  (2.2.17) 

with 

A  =  p  X  q  matrix  (p  >  q  since  oveidetermined) 

X  =  q  X  1  vector  (in  Prony's  method,  these  are  the  aq's) 
i2  =  p  X  1  vector. 

Next,  the  normal  equations  for  the  complex  matrix  A  are  expressed  as: 


AHAx  =  AHj2,  (2.2.18) 

Substituting  QR  for  A  in  this  system  yields  the  following: 

(R"Q")QRx  =RHQ»ii  (2.2.  J  9) 

(Q'lQ)Rx=Q”b  (2.2.20) 

Rx  =  Q'’b,  (2.2.21) 


So, 


X  =  R-IQH  b  (2.2.22) 

The  advantage  to  using  this  approach  over  the  normal  equations  is  evident  when  the 
columns  of  A  are  barely  uncorrelated.  The  ensuing  calculation  of  A^A  will  amplify  round 
off  errors  due  to  this  matrix  being  ill  conditioned.  The  orthogonal  matrix  decomposition 
approach  avoids  this  problem  and  the  round  off  error  accumulations  are  at  a  tninimum[41]. 

Specification  of  a  model  order  is  inherent  in  the  first  step.  The  determination  of  a 
"good"  order  is  complicated  by  using  an  all  pole  (also  known  as  an  autoregressinve  [AR]) 
model  to  represent  a  pole  zero  (also  knowns  as  ARMA)  process.  Even  if  the  number  of 
poles  are  known  apriori,  use  of  the  exact  number  of  system  poles  may  not  yield  a  good 
result .  If  the  model  order  is  underdetermined,  the  "spectrum"  will  be  smooth  and 
smeared.  If  the  order  is  overdetermined,  the  model  is  likely  to  have  spurious  peaks[32]. 
The  empirical  rule  used  by  the  signal  processing  community  of  overestimating  the  model 
order  is  supported  by  two  assumptions.  First,  underestimation  of  model  order  will  not 
identify  true  poles  while  overestimation  will  tend  to  identify  these  poles.  Although  the 
model  is  forced  to  find  parameters  to  fit  the  specified  system,  the  energy  of  the  arbitrary 
poles  (ie,  those  in  excess  of  actual  system  order)  is  quite  small.  Second,  the  presence  of 
noise  in  the  data  may  be  modelled  as  zeros[31].  An  ARMA  process  represented  by  an  AR 
model  requires  overspecification  of  model  order  (actually,  the  bias  of  the  estimation 
decreases  as  the  order  increases). 

Various  researchers  have  suggested  analytical  methods  to  estimate  model 
order[33,42-44].  Akaike[42]  has  suggested  a  final  prediction  error  method  and  a  cost 
minimization  method  in  which  a  cost  is  assigned  for  extra  coefficients  which  do  not  reduce 
model  order.  Criticism  of  the  final  prediction  method  is  that  it  yields  too  low  a  model  order 
while  the  cost  minimization  method  is  said  to  have  statistical  inconsistencies.  To  date,  there 
is  no  common  approach  for  identifying  a  good  model  order.  In  the  algorithm  used  licre, 
the  model  order  specification  is  left  to  the  user.  A  first  approach  used  an  singular  value 
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decomposition(SVD)  of  the  signal  matrix  to  obtain  singular  values  of  the  system.  A 
sudden  decrease  in  the  magnitude  of  the  singular  values  was  the  breakpoint  of  estimating 
the  actual  modes  of  the  system.  Two  problems  with  this  approach  are  that  the  "breakpoint" 
is  not  clearly  defined  as  noise  is  introduced  and  the  optimum  model  order  is  still  not 
specified  by  examining  the  singular  values.  A  more  subjective  approach  made  possible  by 
the  interactive  nattire  of  the  algorithm  was  to  increase  the  order  and  examine  the  outputs. 
When  the  number  and  values  of  the  wavenumbers  of  propagating  modes  (identified  by  low 
damping  and  high  amplitude)  stopped  changing  as  order  increased,  the  model  was  said  to 
have  sufficient  order. 

Polynomial  rooting  is  the  second  step  of  the  Prony  process.  Since  this  is  done 
numerically,  a  robust  algorithm  must  be  used.  The  nonlinear  aspect  of  the  parameter  fitting 
is  located  in  this  portion  of  the  algorithm;  slight  errors  in  the  coefficients  may  result  in 
significant  changes  in  the  roots[41].  For  example,  the  polynomial  x'^  -  lOx^  +  35x2 . 50^ 

+  24  has  roots  of  (x  -  l)(x  -  2)(x  -  3)(x  -  4).  A  change  of  0.5%  in  the  second  coefficient 
yields  a  polynomial  of  x4  -  lO.OSx^  +  35x2-  50x  +  24,  which  has  roots  of  (x  -  0.992)(x  - 
2.340  -  j0.2269)(x  -2.340  +  j0.2269)(x  -  4,378).  The  choice  of  a  complex  rooting  routine 
by  Jenkins  and  Traub  and  double  precision  calculations  are  the  tools  used  to  reduce  errors 
in  this  section  of  the  algorithm[45].  While  the  three  stage  rooting  algorithm  has  performed 
well,  the  sensitivity  of  the  roots  of  a  polynomial  to  the  coefficient  values  indicates  the 
choice  of  algorithms  in  step  one  may  be  the  major  contributing  factor  to  the  accuracy  of  the 
final  model. 

The  third  and  final  step  of  the  routine  is  the  solution  of  llie  overdetennined  version 
of  the  linear  system  expressed  in  equation  (2.2.3).  llie  minimization  of  least  scpiarc  error 
technique  is  used;  tlie  result  is  the  normal  equations  approach  outlined  in  equation  (2.2.1 5). 
A  QR  decomposition  or  the  usual  least  squares  teclinique  (x  -  (A'^A)'k\l*  b  )  may  be 
u.sed(46).  The  algorithm  exercised  in  the  thesis  research  used  the  second  metlK'd  and 
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incoiporated  a  Cholesky  decomposition  scheme  to  exploit  the  Hermitian  symmetry  of 

aha. 
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2.3  Application  of  Prony*s  Method  to  the  Shallow 
Water  Waveguide 

In  the  section  on  underwater  acoustic  theory,  an  expression  for  the  far  field  acoustic 
pressure  as  a  sum  of  normal  modes  was  developed.  While  the  algorithm  for  using  Prony's 
method  to  fit  parameters  to  a  deterministic  model  was  outlined  in  the  last  section,  the  task 
of  transforming  the  pressure  field  to  a  suitable  form  remains.  We  will  now  develop  a 
model  of  the  shallow  water  waveguide  in  a  format  which  fits  that  assumed  in  the  last 
section. 

From  equation  (2.1.8),the  pressure  field  may  be  expressed  as[13]: 

P 

p(r,z)  =  jtc  u*(zo)  Uq(z)  (kqr)  +  I(r)  (2.3.1) 

q=l 


The  far  field  contribution  of  the  continuum,  I(r),  may  be  neglected  and  the  asymptotic 
approximation  for  the  Hankel  function  of 

H^?(kq*-)  =  ^  ■  P  (2.3.2) 


is  substituted  into  equation  (2.3.1)  to  yield: 


p 

p(r,z)  =  jic  ^^^aq  u*(zo)  uqfz)  eJ 


(kqr-  j)  +  aqr  (2.3.3) 


q=i 

Since  the  pressure  field  is  measured  on  a  horizontal  array,  p(r,z)  will  be  expressed  as  p(r) 
and  the  depth  dependence  will  be  incorporated  into  a  constant,  Aq,  in  the  following  manner: 


P 

p(r)  =  "^expKaq  +  jkq)r  +  j0q]  (2.3.4) 

«f=l 


The  data  available  for  processing  are  actually  samples  of  the  pressure  field  rather 


than  the  continuous  pressure  field  itself  The  di.screte  samples  allow  r  in  equation  (2.3.4) 


to  be  replaced  by  nT  in  which  T  is  the  sampling  range.  This  assumes  equally  spaced  data 
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points;  in  practice,  a  cubic  spline  interpolation  scheme  is  used  to  ensure  proper  spacing  of 
the  input  data.to  the  algoriihm.  Experimental  results  indicate  the  sampling  must  be  dense 
enough  to  meet  the  Nyquist  criteria  to  avoid  aliasing.[25]. 

The  last  step  in  modeling  the  pressure  field  in  a  Prony  format  is  the  elimination  of 

the  i  term  in  equation  (2.3.4).  The  data  is  multiplied  by  the  Vr  to  yield  the  model: 

Vr 

y[n]  =  p(r)Vr  =  ^  Aq  exp[(aq  +  jk:q)nT  +  jGq]  (2.3.5) 

q=0 

where 

* 

Aq  =  Ijtc  aq  Uq(zo)  Uq(z)l  =  amplitude 

oq  =  Im(kq)  =  damping 

kq  =  Re(kq)  =  modal  eigenvalues 
7C  * 

0q  =  ^  +  Z  (jjc  aq  Uq(zo)  Uq(z)]  =  initial  phase 


This  model  may  be  further  compacted  in  the  form  of  equation  (2.2.2) 

y[n]  =  ^vqz^  (2.3.6) 

q=0 


where 


Vq  =  Aqexp(j0q) 

Zq  =  exp(((Xq  +jkq)T] 

Once  the  Vq  and  Zq's  are  obtained  through  Prony’s  method,  the  parameters  of  the  model  are 
generated  via: 


tan'i 


(2.3.7) 


aq  = 


In  Izfjl 
T 


Aq  —  I  v'q  I 
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An  examination  of  the  output  of  the  program  for  a  simple  model  is  useful  as  a  guide 
to  gaining  physical  insight  into  the  analysis  results.  The  algorithm  described  in  section  2.1 
was  incorporated  in  a  FORTRAN  program  called  PRAWNS  (PRony  Analysis  of 
Waveguide  for  Nominal  Spectrum).  A  sound  speed  profile  (summarized  in  figure  2.3.1) 
was  developed  and  the  corresponding  complex  pressure  field  versus  range  data  was 
generated  by  SNAP[47].  The  sound  speed  profile  is  a  simplified  version  of  the 
experimentally  determined  profile  of  the  water  column  and  sediment  layers  off  Nantucket 
Island  in  Massachusetts.  SNAP,  a  normal  mode  acoustic  propagation  modelling  program, 
generates  a  far  field  approximation  of  a  pressure  field  from  a  user  defined  profile  and 
frequency  (in  this  case,  220  Hz).  In  addition  to  a  pressure  vs.  range  output,  SNAP 
provides  a  list  of  the  modal  eigenvalues  (the  wavenumbers),  the  normalized  eigenfunctions 
(the  Uq(z)'s)  and  the  attenuation  coefficients.  This  "ground  truth"  permits  a  reasonable 
method  of  examining  the  PRAWNS  output. 

Table  2.3.1  provides  the  SNAP  attenuation  and  wavenumber  outputs  for  this 
profile.  Table  2.3.2  is  a  partial  PRAWNS  output  of  two  different  range  intervals.  (The 
specific  relationship  among  the  model  specifications  listed  in  table  2.3.2  is  covered  in 
Chapter  3).  A  brief  examination  of  the  PRAWNS  outputs  in  table  2.3.2  will  highlight 
some  program  characteristics;  the  information  provided  by  this  simple  analysis  is  rich. 
First,  the  PRAWNS  analysis  provides  a  means  of  identification  of  the  propagating  modes 
in  the  overspecified  Prony  model.  As  expected,  the  propagating  modes  are  marked  by 
high  amplitude  and  low  damping.  As  the  model  order  is  increased,  the  damping  and 
amplitude  remained  constant  for  a  given  range  block.  When  the  order  is  vastly 
overspecified,  there  is  a  small  change  in  Uie  propagating  mode  amplitudes  as  the  energy  is 
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forced  to  be  distributed  among  the  user  specified  order.  In  viewing  the  tabular  data,  high 
amplitude  and  low  damping  proved  to  be  a  good  indicator  as  to  the  actual  order  of  the 
system  (as  mentioned  in  section  2.2.1,  an  analytic  approach  to  obtain  actual  system  order 
may  be  performed  by  an  SVD  of  the  signal  matrix).  Second,  as  mentioned  in  the 
introduction,  the  primary  objective  of  the  analysis  is  to  obtain  accurate  wavenumber  data. 
Comparison  of  the  first  three  indices  of  table  2.3.2  with  the  modes  of  table  2.3. 1  show  a 
wavenumber  match  to  five  significant  digits.  This  is  most  satisfying  given  that  the  aperture 
size  is  90  meters  (30  points  of  data  used).  Third,  the  damping  factors  output  by  PRAWNS 
arc  quite  reasonable  given  the  small  aperture.  It  is  expected  that  the  small  decay  in  energy 
of  a  propagating  mode  over  the  90  meter  range  aperture  would  lead  to  difficulties  in 
accurate  estimation  of  the  attenuation  factor.  There  are  two  ways  to  attempt  to  improve  this 
accuracy.  The  first  method  is  to  increase  the  range  aperture  and  allow  the  energy  of  a  given 
mode  a  greater  distance  to  decay.  The  second  method  uses  the  amplitudes  of  a  given  mode 
as  outlined  in  figure  2.3.2,  The  origin  is  set  to  some  reference  range  and  the  PRAWNS 
algorithm  is  used  on  the  range  interval  R  to  R'.  Since  attenuation  will  cause  an  exponential 
decay  of  the  amplitude  of  a  given  mode  with  respect  to  range,  the  amplitude  of  a  particular 
mode,  A(ro),  estimated  by  PRAWNS  at  a  local  reference  range,  ro,  is  given  by: 

A(ro)  =  A(0)exp[a  rq]  =  A(0)exp[a  RJexpIa  5r]  (2.3.5) 

where 

R»5r 


r  =  0  r  =  R  r  =  R' 

Fig  2.3,2  Amplitude  decay  vs  range 
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While  the  PRAWNS  aperture  may  be  too  small  to  yield  accurate  attenuation  estimates,  the 
use  of  amplitude  estimates  from  seperate  range  blocks  may  permit  a  better  attenuation 
value.  For  a  numerical  example,  consider  the  amplitude  of  mode  1  in  the  1002  - 1 149  m 


and  the  1350  - 1497  m  range  blocks.  Using  a  range  difference  of  1350  -  1002  =  348  m  = 
Ar,  set  up  a  simple  ratio  which  takes  into  account  the  effect  of  attenuation  in  the  model: 


e  a  Ar  ^here  Apq  =  mode  p  amplitude  in  range  block  q.  (2.3.6) 


Substituting, 


1  ,  r0.1044188^ 
““  348  ^\0.1076169j"' 


=  -8.67  E-05 


(2.3.7) 


Comparing  the  results  in  equation  (2.3.7)  to  mode  1  of  Table  2.3.1,  we  see  this  method  is 
not  a  good  choice  in  this  trial.  This  may  be  due  to  poor  amplitude  estimates  or  a  short  range 


interval. 


Table  2.3.1  SNAP  Output  for  Nantucket  Profile 
Mode  No.  Wavenumber  (m- 1 )  Attenuation 

1  0.9074474  7.995E-05 

2  0.8546543  2.422E-04 

3  _ 0.7758756  1.067E-03 

Table  2.3.2  PRAWNS  Output  for  Nantucket  Profile 


Total  No.  of  points:  166  Avg.  block:  50  pts. 

Processing  block:  30  Model  Order:  10 

Overlap:  0  pts.  Samp.  Range:  3.0000  m. 

Starting  Range:  1002.0000  m.  Final  Range:  1 149.0000  m. 


INDEX 

WAVENUMBER 

DAMPING 

AMPLITUDE 

PHASE(RAD) 

1 

0.9074478 

-0.0000796 

0.1076169 

-0.4942496 

2 

0.8546547 

-0.0002414 

0.1850737 

1.5719849 

3 

0.7758779 

-0.0010739 

0.0055408 

-1.1877416 

4 

-0.8612481 

-0.0390015 

0.0000000 

1.5355707 

5 

-0.6309471 

-0.0971785 

0.0000136 

0.4525263 

6 

0.4507092 

-0.1030366 

0.0000168 

- 1 .8049235 

7 

-0.4826431 

-0.1153568 

0.0000122 

0.2122693 

8 

-0.3845713 

-0.1321304 

0.0000194 

-0.9443646 

9 

-0.1642203 

-0.1529004 

0.0000423 

0.2692317 

10 

0.1153333 

-0.1621337 

0.0000536 

2.5536927 

Total  No.  of  points: 

166 

Avg.  block:  50  pts. 

Processing  block:  30 

Model  Order:  10 

Overlap:  0  pts. 

Samp.  Range:  3.0000  m. 

Starting  Range:  1350.0000  m. 

Final  Range:  1497.0000  m. 

INDEX 

WAVENUMBER 

DAMPING 

AMPLITUDE 

PUASE(RAD) 

1 

0.9074471 

-0.0000800 

0.1044188 

-3.0543090 

2 

0.8546540 

-0.0002421 

0.1688871 

-2.0942594 

3 

0.7758755 

-0.0010671 

0.0036974 

-0.0749611 

4 

0.2020056 

-0.0030457 

0.0000000 

-2.3543923 

5 

0.3843271 

-0.0272609 

0.0000000 

0.6840160 

6 

-0.5876253 

-0.0380173 

0.0000000 

2.2279577 

7 

-0.8681644 

-0.0567308 

0.0000000 

-0.643.3010 

8 

-0.0685432 

-0.0614808 

0.0000000 

0.58104X4 

9 

-0.3724043 

-0,0977833 

0.0000001 

0.9746809 

10 

-0.432375 1 

-0.1299475 

0.0000001 

2.0567591 
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As  the  range  between  the  analysis  blocks  increases,  the  decay  in  energy  will  be  greater  and 
the  attenuation  estimate  should  improve.  This  alternate  method  of  estimating  attenuation  is 
also  prone  to  noise  effects  manifested  in  amplitude  estimations. 

Because  of  the  isovelocity  water  coluntn  of  the  Nantucket  profile  of  figure  2.3.2, 
the  mode  shapes  in  the  water  column  are  not  distorted  by  gradients.  The  mode  shapes 
provided  as  part  of  the  SNAP  output  may  be  used  as  a  comparison  to  the  Prony  analysis 
for  a  vertical  array.  This  was  synthesized  by  using  the  sound  speed  profile  of  figure  2.3.1 
but  changing  the  receiver  depth.  The  receiver  depth  was  changed  in  half  meter  increments 
from  0.5  to  13.5  meters.  A  range  aperture  of  375  meters  was  specified  and  the  PRAWNS 
program  was  used  to  analyze  the  modal  structure  at  500, 1000  and  1400  meters  (see  figure 
2.3.3). 


Receiver  depth  changed  in 
0.5  m  depths  from  0.5  to 


13.5  m. 


Fig  2.3.3  Vertical  Array  of  Nantucket  Profile 
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SNAP  and  PRAWNS  mode  shapes  (500  m) 


SNAP  and  PRAWNS  mode  shapes  (1000m) 


SNAP  MODE  1 
SNAP  MODE2 
SNAP  MODES 
Modcl(500m) 
Mode2(500m) 
Mode3(500m) 


SNAP  MODE  I 
SNAP  MODE2 
SNAP  MODES 
Modcl(lOOCIm) 
Modc2(IOOOm) 
ModeS{I000m) 


SNAP  MODE  1 
SNAP  MODES 
SNAP  MODES 
Mode3(  1400m) 
Modc3(l400m) 
Mode3(  1400m) 


Figure  2.3.4  Vertical  Mode  Shapes  for  PRAWNS  and  SNAP 


The  resulting  amplitudes  for  each  mode  were  scaled  to  match  the  normalized  SNAP  mode 
shapes  at  a  particular  depth  (10  meters).  The  resulting  close  agreement  is  evident  in  figure 
2.3.4.  The  larger  aperture  presents  another  opportunity  to  try  the  attenuation  estimation 
method  in  equation  (2.3.6).  With  a  range  interval  of  249  meters  and  receiver  depth  of  10 
meters,  the  results  are  summarized  in  table  2.3.3: 


Table  2.3.3  Attenuation  Estimation  using  Nantucket  Profile 


Mode  No. 

Amplitude  1 

Amplitude  2  Attenuation 

1 

0.2270556 

0.2225807 

-7.994E-05 

2 

0.0687738 

0.0647491 

-2.422E-04 

3 

0.0054065 

0.0041446 

-1.068E-03 

The  results  are  much  closer  to  the  SNAP  results  in  Table  2.3.1. 
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2.4  Summary 

The  cursory  review  of  normal  mode  theory  of  section  2. 1  was  meant  to  emphasize 
key  aspects  of  the  shallow  water  waveguide  problem.  A  general  form  for  the  pressure  field 
was  developed  in  the  form  of  a  modal  sum  and  a  continuum.  A  hard  bottom  example  was 
used  to  demonstrate  the  form  of  the  pressure  field,  the  asymptotic  form  of  the  Hankel 
function  and  the  effect  of  source  receiver  geometry.  The  Pekeris  example  demonrtrated  the 
discrete  spectrum  and  continuous  spectrum  regions  of  more  complex  bottom  models.  The 
distinction  between  the  discrete  and  continuous  spectrum  may  be  set  by  the  wavenumber  of 
a  "basement"  isovelocity  halfspace.  Attenuation  effects  were  incorporated  by  a  small 
imaginary  term  in  the  horizontal  waver  umber,  kn- 

In  section  2.2,  the  basic  three  steps  of  Prony's  method  were  developed.  Upon 
closer  examination  of  these  steps,  some  of  the  options  and  rationale  for  choosing  them 
were  explained.  The  algorithm  used  in  the  remainder  of  the  paper  was  presented;  the 
particular  method  chosen  for  funher  investigation  is  different  from  the  forms  commonly 
found  in  the  signal  processing  literature. 

Section  2.3  tailored  the  modal  sum  representation  of  the  pressure  field  into  a  form 
suitable  for  processing  via  Prony's  method.  Two  applications  of  the  method  were 
presented.  The  first  was  a  straightforward  horizontal  a'my  application  with  a  range 
aperture  of  90  meters  on  data  generated  synthetically  by  SNAP.  This  application 
demonstrated  the  good  agreement  of  wavenumbers  between  PRAWNS  and  SNAP;  in 
addition,  the  distinguishing  characteristics  of  propagating  modes  were  presented.  ITic 
second  application  simulated  a  vertical  array  by  analysis  of  the  same  range  blocks  at 
different  receiver  de;.  ihs.  This  example  demonstrated  use  of  the  PRAWNS  mode 
amplitudes  to  obtain  vertical  eigenfunction  shapes  and  an  aUemalc  method  for  ol.  aining 
mode  attenuation. 
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Chapter  3 

Prony's  Method  Evaluation 


3.1  Depth  Dependent  Green*s  Function 

In  the  previous  chapter,  the  pressure  field  was  represented  by  a  modal  sum  and 
continuum  contribution.  While  the  modal  nature  of  sound  propagation  may  be  observed  in  the 
pressure  field  as  an  interference  pattern,  this  is  not  the  optimum  domain  for  modal  analysis.  As 
the  number  of  modes  increases,  the  interference  pattern  becomes  increasingly  complex.  An 
alternate  descripdon  of  the  modal  information  may  be  obtained  through  a  Hankel  transform  of 
the  pressure  field[48-50].  The  result  of  this  transformation  is  known  as  a  depth  dependent 
Green's  function.  This  section  outlines  the  form  of  the  depth  dependent  Green's  function. 

A  zero  order  Hankel  transform  is  defined  as[21]; 

oo 

H{  F(r)  )  =  f(kr)  =  lF(r)  Jo(kr  r)  r  dr  (3.1.1) 

0 

With  an  inverse  of: 


H{  f(kr)  )  =  F(r)  =  |f(kr)  Jo(kr  r)  k^  dkr 
0 


(3.1.2) 


That  is,  the  Hankel  transform  is  its  own  inverse.  From  equation  (2.16),  recall, 
1  3  3p  92p  -  5(r)  _ 

-—(r •:^)  +  ~  +  k2(z)  p  =  -2— 5(z-zo)  (3.1.3) 

‘  OT  dr  dz^  ‘ 


By  taking  the  Hankel  transform  of  both  sides,  the  equation  yields  an  ordinary  differential 
equation: 


g(kr,z,zo)  =  -  2  6(z-zo) 


(3.1.4) 


In  obtaining  this  equation,  the  following  relationship  is  used: 


f(kr) 


(3.1..S) 


The  solution  to  equation  (3.1.4),  g(kr,z,zo),  is  known  the  depth  dependent  Green's  function. 
An  analytical  solution  to  the  equation  may  be  found  by  application  of  Sturm  Liouville  theory. 
For  a  waveguide  with  the  surface  at  z=0  and  bottom  at  z=h,  the  impedance  boundary 
conditions  may  be  expressed  in  the  form: 

Ai<})'(0)  +  Bi<|)(0)=0  (3.1.6) 

A2(|)'(h)  +  B2  <t)(h)=0 

I  t 

Defining  the  Wronskiain,  W,  as  W  =  ({i^  -  <5>l)  <|)^,  permits  the  solution  to  equation  (3.1.4)  to  be 

expressed  as: 


^  “  {w(zo)}  ‘I’sO^r.z)  <|)b(kr,zo) 

0  <  z  <  zp 

(3.1.7) 

^  ~  {w(zo)}  <l>b(kr.z) 

zp  <  z  <  h 

The  determination  of  the  Green’s  function  now  requires  identifying  and  which 


are  two  solutions  to  a  Sturm  Liouville  problem.  Both  solutions  satisfy  the  homogeneous 
equation  (3.1.4)  and  4)b  satisfies  the  bottom  boundary  condition  while  (])s  satisfies  the  surface 
boundary  condition.  In  terms  of  reflection  coefficients,  for  an  isovelocity  waveguide,  the 
solutions  at  the  surface,  <{)s,  and  bottom,  <{)5,  the  superposition  of  up  and  downgoing  plane 
waves  are[9,26]: 

(t)s  =  A  [e- JYz  +  Rs(kr)  eJYz]  and  (3.1.8) 

(j)b  =  B  [eiY^  +  Rb(kr)  eJY  (2h  -  z)  ] 

Substitution  of  these  into  the  equation  for  the  Green's  function,  equation  (3. 1 .4),  yields: 

j[e- jyiz-zpi  +  Rs  ejY(z^V  +  Rb  cJy^h  (e-iriz+z^)  +  g- JY ] 

^  y[i  -RbRsej2Yh] 

In  the  simple  waveguide  examples  of  chapter  2,  the  surface  reflection  coefficient  is  =  - 1 ;  ie. 
the  upper  boundary  is  treated  as  a  pressure  release.  Tlic  expression  above  Iiolds  equally  well  in 
multilayer  waveguides  in  which  the  upper  layer  boundary  condition  is  not  constrained  to  be  a 


-40- 


pressure  release  meclkanism.  By  the  definition  of  the  Hankel  transform  pair  given  above,  the 
relationship  between  the  depth  dependent  Green’s  function  and  the  pressure  field  is  then[14]: 


p(r,z,zo)  =  J  g(k;r,z,zo)  JoCkrO  kr  dkf  and  (3.1.10) 

0 


g(kr,z,zo)  =  /  p(r,z,zo)  Jo(krr)  r  dr 
0 


The  transformation  of  the  pressure  field  to  a  depth  dependent  Green's  function  is  a 
particularly  useful  tool  in  examining  the  modal  behavior  of  the  waveguide.  The  depth 
dependent  Green's  function  (which  will  be  called  the  Green’s  function  hereafter),  is  a  function 
of  the  horizontal  wavenumber,  kr,  and  the  source  and  receiver  depths,  zq  and  z  respectively.  It 
is  in  the  examination  of  the  horizontal  wavenumber  spectra  that  the  influence  of  the  boundaries 
on  the  modal  structure  is  easily  observed. 

In  chapter  2,  the  pressure  field  was  also  expressed  as  a  sum  of  a  modal  portion  and  a 
continuous  contribution.  Current  studies  of  the  normal  mode  approach  have  shown  this 
decomposition  carries  over  to  the  Green’s  function  in  the  following  manner[25]: 


N 


(1) 

P(r)  =  j  tt  X  a n  k  n,  Hq  (k^r)  +  p  (r) 


n  =  1 


I  ,  I 

8(k,)=  X 


n=l  ''  r 


Oa  If 

+g^(k^) 


(3.1.11) 


(3.1.12) 


That  is,  the  modal  portion  of  the  pressure  field  is  directly  related  to  the  modal  portion  of  the 
Green's  function  through  a  zero  order  Hankel  transform.  The  coefficient,  an,  is  the  residue 
(from  the  Cauchy  residue  theorem)  for  a  given  kr;  that  is,  the  coefficient  is  defined  as  an  = 
(kr-km)  g(kr).  An  alternate  expression  for  the  coefficient  in  terms  of  the  vertical 
eigenfunctions,  Un,  isf25]: 
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3n(kr)  - 


(3.1.13) 


Un(zo)  Un(z) 

P(Z0)  km 


For  a  Pekeris  waveguide,  the  specific  coefficients  may  be  found  using: 

sin(Yoz)  sin(Yozo) _ 


Sn(kr)  -■ 


,  ,  ,  sin(2Yoh)  .  b  sin2(Yoh) 
krl  h - ; - +  j 


4Y0 


2yi 


) 


kr  =  km 


(3.1.14) 


where 


b  = 


YO 

j  Yl  h 


kr  =  km  the  subscripts  0  and  1  denote  the  water  column  and  half 


space  respectively. 

There  are  two  additional  items  which  should  be  noted  with  regard  to  the  Green's 

function.  The  first  characteristic  concerns  the  rational  form  expression  for  the  Green's 

function.  In  equation  (3.1.12),  the  modal  portion  of  the  Green's  function,  gm,  is  given  as; 

N 

gm=  X  0^"^%  ■  (3.1.15) 


I 


(^r  - 


n=l 


where  both  the  numerator  and  denominator  of  the  expression  have  zeroes  at  locations  other 
than  zero  and  infini^.  In  signal  processing  parlance,  this  is  a  pole- zero  or  ARMA 
(autoregressive  moving  average)  model.  The  zeroes  of  the  denominator  correspond  to 
singularities  and  are  called  the  poles  of  the  system.  In  the  above  expression,  the  poles  are 
expressed  in  terms  of  the  horizontal  wavenumber  and  are  located  at  kr  =  ±  km(as  expected 
since  the  Green's  function  is  an  even  function).  The  zeroes  of  the  numerator  are  the  nulls  of 
the  system  and  are  found  through  the  identification  of  the  coefficients  an(kr),  which  may  be 
done  analytically  or  numerically. 

The  second  fact  regarding  the  Green's  function  is  based  on  the  physical  interpretation 
of  equation  (3.1.9).  The  poles  of  the  equation  are  functions  of  the  waveguide  environment 
since  the  numerator  is  a  function  of  Rb,  Rs.  h  and  kr  (recall  the  horizontal  wavenumber,  kp,  is 
related  to  the  vertical  wavenumber,  y,  by  kr  =  yk^  -  y^  )•  3'he  numerator  of  the  Green's 
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function  has  zeroes  which  are  determined  by  the  waveguide  environment  gnd  the  geometry  of 
source  and  receiver.  As  the  source  and  receiver  positions  are  varied  within  a  given  waveguide, 
we  expect  a  shift  in  the  zero  locations  and  a  corresponding  change  in  the  null  locations  on  a 
spectral  plot  of  the  Green's  function.  At  certain  source-receiver  geometries,  the  zero  will 
cancel  a  pole  and  the  mode  will  not  propagate. 
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3.2  Prony  Energy  Spectral  Density 

In  the  last  section,  the  modal  behavior  of  the  waveguide  was  available  for  examination 
through  the  Green's  function.  The  Green's  function  is  related  to  the  pressure  field  through  a 
zero  order  Hankel  transform.  Prony's  method  is  a  model  parameter  estimation  approach  rather 
than  a  spectral  estimation  technique.  Since  Prony's  method  estimates  the  system  poles  ( the 
km)  as  one  of  the  model  variables,  we  wish  to  develop  a  spectral  plot  which  incorporates  this 
and  the  other  parameter  values. 

The  first  step  in  defining  a  spectrum  for  the  method  is  an  assumption  of  the  behavior  of 
the  data  outside  the  processing  interval.  Marple  presents  three  possible  spectra  which  start 
from  different  interpretations  of  data  characteristics[31].  In  this  study,  the  data  is  assumed  to 
have  even  pressure  characteristics  with  respect  to  the  origin,  that  is,  p(r)  =  p(-r).  Physically, 
this  may  be  justified  by  considering  the  symmetry  assumptions  made  in  the  analytical  solution 
to  the  waveguide  problem.  In  chapter  2,  a  circular  symmetry  assumption  was  invoked  to  allow 
removal  of  0  dependency  through  integration.  In  terms  of  measurements,  if  we  place  our 
spatial  array  at  a  horizontal  range  xq  from  the  origin,  we  would  expect  the  same  measurements 
if  the  array  was  located  at  -xq.  Simply  put,  this  assumption  requires  the  effects  of  attenuation, 
spreading,  etc.  to  be  the  same  on  either  side  of  the  origin. 

This  symmetry  is  incorporated  in  the  spectrum  by  assuming  a  two  sided  function  which 
is  defined  as: 


The  z-transfomi  of  this  model  is[51]: 
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Using  the  definitions  of  the  roots  Zq  and  z*  : 

Zq  =exp(aqT+j  kqT) 
(z*)-l=exp(-aqT  +  jkqT) 


(3.2.2) 


(3.2.3) 

(3.2.4) 


and  assuming  the  spectral  radius  of  the  poles  is  less  than  one,  a  discrete  time  Fourier  transform 
(DTFT)  may  be  found.  The  DTFT  is: 

Si(kr)  =TY(z)L  =  exp(jkrT)  (3.2.5) 

P-1 

_ T  [exp  (OqT)  -  exp(-aqT)  ]exp(j[krq  -  kr]T) _ 

1-  [exp  (OqT)  -  exp(-aqT)]exp(ji;krq  -  kJT)  +  exp(j2[krq  -  kr]T) 
q-0 


The  energy  spectral  density  (ESD)  is  found  by  the  magnitude  squared  of  the  DTFT,  ie  S(kr)  = 
ISi(kr)|2.  We  will  refer  to  this  function  as  the  ESD  in  the  rest  of  the  thesis. 

The  z  transform  of  the  two  sided  function  demonstrates  a  characteristic  of  the  Prony 
method;  the  model  is  an  all  pole  system.  Rewriting  the  z  transform  of  equation  (3.2.1)  in 


terms  of  a  common  denominator: 

P-1 


Y(z)  = 


1  -  (ZqZ)-l  -  (1  -  ZqZ-1) 

(1  -  ZqZ-1)  (1  -  (z*z)-l) 


(3.2.6) 


[^q  - 

(1  -  ZqZ-1)  (1  -  (z*z)-l) 


(3.2.7) 


In  the  above  representation,  it  is  readily  apparent  the  numerator  of  the  z  transform  has  no 
zeroes  except  at  zero  and  infinity  for  all  q. 
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The  all  pole  nature  of  the  Prony  ESD  dictates  the  appearance  of  the  spectral  plot.  The 
singularities  of  the  system,  the  poles,  yield  a  typical  spectrum  which  has  sharp  peaks.  In 
contrast,  an  all  zero  (MA  or  FIR)  filter  has  sharp  nulls.  The  pole-zero  model  (ARMA)  has  both 
sharp  nulls  and  peaks  in  the  spectral  plot[51]. 


Fig  3.2.1  Comparison  of  typical  AR,  MA  and  ARMA  spectra 
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In  the  last  section,  the  pole-zero  nature  of  the  Green's  function  was  evident.  Indeed,  a 
plot  of  the  Green's  function  for  synthetically  generated  data  confirms  this  characteristic  (see, 
for  example,  fig  3.3.3).  Instead  of  taking  a  Hankel  transform  of  the  pressure  field,  the  Prony 
BSD  is  achieved  by  a  DTFT  of  the  square  root  of  the  range  times  the  modal  portion  of  the 
discrete  pressure  field.  While  the  ESD  of  the  resulting  all  pole  model  not  the  Green's  function, 
it  is  related  to  it  The  Wold  decomposition  theory  states  that  we  may  model  an  ARMA  or  MA 
filter  with  an  infinite  order  AR  filter[32].  Given  the  constraint  of  a  finite  filter,  the  relationship 
between  the  AR  filter  coefficients  and  the  ARMA  coefficients  and  the  AR  filter  length  should  be 
explored.  Consider  the  Green's  function  as  an  ARMA  model  of  s  poles  and  r  zeroes.  Then, 
using  the  analogy  of  a  filter  driven  by  a  unit  impulse; 


H(z)  = 


2b[m] 
B(z)  m=0 
D(z) 


,-m 


(3.2.8) 


2a[ql  z-q 

q=0 


with  b[0]  =  a[0]  =  1.  To  match  this  with  an  ail  pole  filter  specified  by: 

^g[w]  z-'^ 
w=0 

requires  H(z)  =  V(z).  Substituting  the  above  expressions: 

^  =>  B(z)  G(z)  =  D(z)  (3.2. 10) 

The  multiplication  in  the  z  domain  is  equivalent  to  convolution  in  the  time  domain,  so: 

b[n]*g[n]=d[n]=>  Z(b[kl  g[n-k])  =  d[n]  (3.2.11) 

k=-oo 


Since  there  are  s  poles  in  the  ARMA  model,  d[n]  =  0  for  n  >  s  and 


g[n]  =  -  Z(b[k]  g[n-k]) +d[n] 
k=l 

=  -  i(b[k]  g[n-k]) 
k=l 


0  <  n  <  s 


n  >s 


(3.2.12) 
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This  expression  assumes  g[n]  is  causal  and  g[0]  =  1.  The  validity  limits  for  the  expressions 
are  due  to  the  number  of  ARMA  poles  while  the  number  of  ARMA  zeroes  defines  the 
convolution  length.  In  the  AR  approximation,  these  g[n],  the  all  pole  filter  coefficients,  are  set 
equal  to  zero  for  n  >  r+s.  As  shown  above,  the  poles  and  zeroes  of  the  ARMA  filter  influence 
the  output's  first  s  points.  After  that,  only  the  poles  influence  the  system  output.  To  determine 
the  ARMA  parameters  given  an  AR  approximation,  the  above  equations  may  be  used  to  find 
the  numerator  and  denominator  terms  using  r+s+1  output  points: 


/  c[s]  c[s-l] 

c[s+l]  c[s] 

\c[s+r-l]  c[p+q-2] 

and 

a[n]  =  c[n]  +  ib[k]c[n-k] 
k=l 


c[s-r+l]  \ 

^[1]\ 

/c[p+l]\ 

c[S'r+2] 

b[2]  1 

c[p+2] 

c[p]  J 

\mJ 

Vc[p+q]/ 

(3.2.13) 


(3.2.14) 


In  other  words,  the  expressions  for  g[n]  are  forcing  these  coefficients  to  match  the  first  r+s+1 
coefficients  of  the  infinite  length  inverse  ARMA  polynomial  X(z)  where 

oo 


X(z)  = 


_1 _ 

H(z)  -  B(z) 


k=0 


(3.2.15) 


The  last  section  contained  a  description  of  the  Green's  function  as  a  superposition  of  a  modal 

and  continuous  portion.  The  modal  sum  in  the  Green’s  function  may  be  expressed  by [25]; 

N 

gm(kr)  =  \7  (3.2.16) 


I 


where 


n=l 

Un(zo)  Un(z) 
an  — - 

p(zo)  krn 


This  may  be  cast  in  the  form  of  an  ARMA  model  with  kp  replaced  by  z.  However,  the 
expression  for  the  coefficients  an  is,  in  general,  not  trivial  and  requires  substitutions  and/or 
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series  expansions  to  yield  a  polynomial  form  in  the  numerator.  A  direct  comparision  of  the 
Green's  function  and  Prony  ESD  is  not  useful;  the  quantities  are  not  the  same  since  the 
Green's  function  is  an  ARMA  process  and  the  Prony  ESD  is  an  AR  process.  The  poles  of  the 
Green's  function  will,  however  mathc  those  of  the  Prony  ESD;  the  ESD  will  also  yield  relative 
energy  levels  of  the  system  modes.  The  Green's  function  is  related  to  the  pressure  field  by  a 
zero  order  Hankel  transform  while  the  Prony  ESD  is  related  to  the  square  root  of  the  range 
times  the  pressure  field  by  a  discrete  time  Fourier  transform. 

The  ESD  uses  all  of  the  information  generated  by  the  Prony  estimation  algorithm.  This 
presentation  has  two  distinct  advantages.  First,  the  tabular  form  is  transformed  to  a  graphical 
representation.  The  total  effect  of  all  model  parameters  is  summarized  in  a  concise  output 
which  allows  easy  assimilation  of  the  algorithm  output.  Second,  the  ESD  is  a  tool  which 
allows  comparision  of  range  blocks.  The  present  algorithm  has  no  constraints  on  continuity  of 
parameter  values  between  subsequent  range  blocks;  each  block  is  evaluated  as  a  "stand  alone" 
entity.  Changes  in  bathymetry  and  waveguide  boundaries  may  be  observed  in  shifts  in  the 
ESD  peaks  and  levels  between  blocks  of  interest.  This  range  dependent  performance  is  funher 
explored  in  the  next  chapter. 
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3.3  Algorithm  Performance  in  Range  Independent 
Waveguides 

In  this  section,  we  will  explore  the  application  of  the  Prony  algorithm  to  a  synthetically 
generated  pressure  field.  The  algorithm  explored  was  outlined  in  chapter  two;  in  brief,  it  is  a 
three  step  process.  The  first  step  consists  of  obtaining  a  least  squares  fit  of  the 
polynomial  coefficients  for  a  specified  model  order.  After  using  a  QR  decomposition  to  obtain 
these  coefficients,  a  numerical  rooting  program  finds  the  roots  of  the  polynomial.  The  last  step 
is  the  LS  fit  of  the  data  to  the  remaining  model  parameters.  The  parameters  specified  by  the 
users  include  the  starting  and  stopping  range,  the  model  order,  the  processing  block  size, 
processing  block  overlap  and  averaging  block  size.  The  relationship  of  the  last  three  variables 
is  outlined  below  and  in  figure  3.3.1. 


The  processing  block  is  the  number  of  points  used  for  each  three  step  iteration  of  the  method. 
The  processing  block  overlap  allows  the  user  to  set  the  starting  and  stopping  ranges  of  each 
block  so  that  adjacent  blocks  may  overlap  one  another.  The  averaging  block  size  pcnnits  the 
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parameter  values  within  a  given  range  interval  to  be  averaged  together.  In  figure  3.3. 1,  an 
averaging  block  is  shown  with  three  processing  blocks  overlapping  by  one  third.  It  should  be 
emphasized  that  this  is  not  an  averaging  of  spectra  as  is  usually  encountered  in  the  signal 
processing  literature.  Rather,  the  averaging  is  among  the  parameter  values;  it  takes  place 
before  the  ESD  is  calculated. 

Before  exploring  the  utility  of  the  ESD  in  evaluating  the  Prony  algorithm  performance, 

we  will  describe  another  analysis  tool  which  may  be  used  to  obtain  a  quantitative  assessment  of 

the  modelling  process.  As  shown  in  figure  3.3.2  below,  the  observed  data  (the  square  root  of 

the  range  times  the  sampled  pressure  field)  is  assumed  to  be  the  output  of  an  ARMA  filter.  The 

all  pole  filter  of  Prony's  method  uses  the  observed  data,  y[n],  to  estimate  the  Prony  filter 

coefficients,  Vfz).  Given  the  all  pole  filter,  V(z),  we  want  to  develop  a  method  of  assessing 

how  well  the  actual  data  fits  the  all  pole  model.  If  s[n],  the  output  of  the  Prony  filter,  was  used 
as  input  of  the  inverse  filter,  G(z)  =  the  output,  r[n],  would  be  a  unit  impulse[51]. 

Instead  of  using  s[n]  as  input  to  G(z),  the  observed  data  is  used.  If  H(z)  was,  in  fact,  all  pole, 

and  matched  by  V(z),  then  the  output  of  G(z)  would  be  an  impulse  of  height  equal  to  the  gain 

term  of  H(z)  (because  V(z)  assumes  a  gain  of  1).  If  H(z)  is  an  ARMA  process  and/or  there  is 

noise  in  the  system,  then  using  y[n]  as  input  to  the  inverse  filter  will  yield  a  data  sequence 

which  has  nonzero  terms  at  other  than  the  origin.  The  residue  is  the  normalized  total  energy  of 

this  sequence:  R  =  X  *‘*[n].  The  smaller  the  residue,  the  better  the  actual  data  fits  an 
n=l 

all  pole  model  since  there  is  less  energy  in  the  residual.  The  last  statement  requires  a 
qualification  for  completeness.  The  residue  information  yields  a  quantitative  measurement  of 
observed  data  fit  to  an  all  pole  model;  it  makes  no  statement  regarding  the  accuracy  of  the 
parameter  values  to  the  actual  waveguide  values. 
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Actual  ARMA  process  All  pole  model 


Fig  3.3.2  Residue  determination 


With  these  caveats,  it  is  acceptable  to  add  the  residue  to  the  performance  evaluation 
tools.  The  four  tools  which  will  be  used  are  the  BSD,  residue,  wavenumber  vs.  range  and  pole 
positions  in  the  complex  z  plane.  Each  of  these  serves  to  highlight  specific  characteristics  of 
the  waveguide  performance.  The  ESD  is  best  used  to  summarize  the  modal  behavior  of  the 
system  including  range  dependent  features  (which  will  be  covered  in  chapter  4).  The  residue 
and  wavenumber  plots  are  useful  for  determining  when  a  "good"  system  order  is  reached  and 
the  pole  plots  assist  in  identifying  actual  modes  from  the  arbitrary  system  poles. 

One  of  the  classical  difficulties  facing  the  researcher  implementing  Prony’s  method  is 
the  specification  of  a  model  order.  This  parameter  is  distinct  from  the  system  order.  The 
system  order  (the  number  of  modes  actually  propagating)  may  be  obtained  by  inspection  or 
numerically.  The  Green's  function  (which  is  generated  by  numerically  Hankel  transforming 
the  pressure  field)  may  be  used  to  estimate  the  number  of  system  modes.  However,  large 
apertures  are  necessary  for  accurate  discrete  Hankel  transform  results.  To  numerically  extract 
the  number  of  system  modes,  an  SVD  of  the  signal  matrix  is  the  recommended 
procedure[36, 37 ,43-44].  Ordering  of  the  singular  values  by  magnitude  will  illustrate  the  drop 
in  singular  value  magnitude  for  singular  values  greater  than  the  actual  system  order  ( ic,  if  the 
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mcxicl  order  is  ten  and  the  system  order  is  six,  the  first  six  largest  singular  values  will  be 
significantly  higher  than  the  remaining  singular  values).  The  addition  of  noise  to  the  system 
complicates  both  of  these  procedures. 

The  specification  of  a  "good"  model  order  is  more  of  an  art.  If  the  order  specified  is 
too  low,  the  resulting  parameter  estimation  will  be  poor  as  the  system  energy  is  constrained  to 
be  distributed  among  the  specified  order.  Empirically,  this  situation  results  in  the  identification 
of  some  of  the  actual  modes,  generally  the  strongest  (which  arc  determined  by  source-receiver 
geometry).  The  sytem  description  is,  nonetheless,  incomplete.  If,  on  the  other  hand,  the 
model  order  is  too  high,  the  output  may  contain  spurious  peaks  and/or  the  signal  matrix  may  be 
singular[32].  While  the  noise  and  non  modal  components  of  the  energy  field  ( the  continuum) 
may  prevent  the  matrix  from  actually  becoming  singular,  the  matrix  becomes  more  ill 
conditioned  as  order  increases.  This  may  make  numerical  determination  of  the  actual  system 
order  easier  since  the  larger  singular  values  get  larger  and  the  smaller  singular  values  get 
smaller  near  singularity[40,41].  Since  we  do  not  specifically  use  the  actual  system  order,  this 
small  positive  feature  of  largely  overspecified  models  is  negated  by  the  problems  associated 
with  decomposition  of  the  ill  conditioned  matrix.  The  spurious  peaks  are  usually  of  low 
energy  and  have  little  effect.  In  some  trials,  the  small  but  finite  energy  in  these  poles  robbed 
the  valid  modes  of  energy  resulting  in  a  poorer  fit  of  data  to  an  all  pole  model.  The 
development  of  an  analytical  method  to  determine  an  appropriate  model  order  is  a  current 
research  effort.  In  a  recent  paper,  Braun  and  Ram  describe  an  SVD  approach  which  yields 
effective  results  esfiecially  in  noisy  situations[6]. 

The  approach  used  in  this  study  to  examine  the  effect  of  model  order  on  parameter 
estimation  was  iterative  in  nature.  Two  different  schemes  were  used.  In  the  first,  the  model 
order  was  increased  and  the  total  residue  for  the  range  interval  was  obtained.  The  residue 
dropped  to  a  plateau  after  a  certain  amount  of  overdetermination.  The  second  scheme  examined 
parameter  "wander"  with  respect  to  increases  in  model  order.  Both  pole  plots  and  wavcnumticr 
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vs.  range  plots  were  examined  with  increases  in  model  order.  Initially,  if  the  model  was 
underspecified  in  order,  the  wavenumbers  and  other  parameters  were  found  to  change  slightly 
as  the  order  was  increased.  At  a  certain  point  in  overspecification,  the  wavenumber  change 
was  negligible  for  increase  in  model  order.  The  pole  plot  proved  to  be  more  effective  in 
determining  actual  modes  rather  than  in  assisting  in  parameter  wander  identification. 

In  exercising  the  algorithm,  two  bottom  models  were  chosen.  The  first  was  a  Nantucket 
type  bottom  with  source  frequencies  of  140  and  220  Hz.  The  second  bottom  profile  was 
similar  to  the  type  of  bottom  found  off  the  coast  of  Corpus  Christi  with  source  frequencies  of 
50  and  140  Hz.  The  bottom  models  were  chosen  because  they  approximated  the  experimental 
environments.  The  sound  fields  were  generated  using  the  SAFARI  code  which  utilizes  a 
propagator  matrix  approach[55,56].  Unlike  SNAP,  continuum  contributions  are  addressed  by 
SAFARI.  Because  of  the  continuum  contribution,  we  expect  the  aU  pole  model  not  to  fit  the 
pressure  field  as  well  as  a  SNAP  generated  data  set.  The  continuum  effects  are  near  field;  at 
long  ranges,  their  effect  should  be  quite  small. 

In  the  Nantucket  model,  the  bottom  profile  of  figure  2.3.1  was  used  with  an  attenuation 
of  0.07  dB/A,.  The  three  propagating  modes  are  well  defined  and  widely  spaced;  this  bottom  is 
not  particularly  sensitive  to  model  parameter  variations  in  model  order,  processing  block, 
averaging  or  overlap.  Nonetheless,  the  bottom  model  serves  to  illustrate  some  basic 
performance  aspects  of  the  algorithm  in  a  "realistic"  (but  still  noisefree)  environment.  The  ESD 
provides  a  useful  tool  for  the  overall  modal  behaivor  of  the  waveguide.  In  this  respect,  the 
ESD  is  comparable  to  the  Green's  function  (figure  3.3.3)  since  it  indicates  energy  in  the 
various  modes.  In  figure  3.3.4,  a  small  aperture  model  is  used  to  estimate  the  parameters  for 
the  ESD.  The  total  number  of  points  needed  for  each  range  block  is  15  points  at  3.2  meter 
spacing. 
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Fig  3.3.4  ESD  of  PRAWNS  output  on  Nantucket  profile  at  4000  m 
(15  points,  0  avg,  3.2  m  spacing,  model  order  7) 

The  accuracy  of  the  PRAWNS  algorithm  output  vs  the  Green's  function  pole  locations 
is  best  seen  in  table  3.3.1.  The  Green's  function  pole  locations  were  found  through  use  of  a 
peak  finding  routine.  The  PRAWNS  wavenumbers  are  the  tabular  values  of  the  model 
specified  in  figure  3.3.4  at  4000  meters.  The  range  of  4000  meters  refers  to  the  first 
processing  block  which  contains  that  user  specified  range;  the  4000  meter  point  may  occur 
anywhere  within  this  processing  block.  This  manner  of  specifying  a  range  will  be  the 
convention  in  the  rest  of  this  thesis. 

Table  3.3.1  Comparision  of  PRAWNS  and  Green's  function  pole  locations 

(Nantucket  profile) 

Mode  Horizontal  Wavenumber  (m'^) 

PRAWNS  Green's  Function 

1  0.9077464  0.907701 

2  0.8569943  0.856988 

3  0.7822649  0.782204 


The  Nantucket  profile  provides  the  opportunity  to  demonstrate  the  utility  of  the  pole 
plots.  This  is  a  Z  plane  representation  of  the  wavcnumlxir  and  damping  parameters  using  the 
infonnation  generated  as  follows: 


-.56 


(3.3.1) 


Pole  magnitude  =  VZreal^  +  Zimag^ 

Pole  angle  =  tan* 

where  Zreal  =  *  T  cos(  wavenumber  *  T) 

2^ag  =  *  T  sin(  wavenumber  *  T) 

for  each  damping  and  wavenumber  term.  Identifying  valid  modes  among  the  arbitrary  poles 
estimated  by  the  overspecified  system  is  performed  by  recognizing  that  the  valid  modes  won't 
change  with  variations  in  order  while  the  arbitrary  modes  will.  Figure  3.3.5  illustrates  the 
model  poles  for  various  orders  and  the  overlap  situation  which  allows  pole  identification. 

The  use  of  residues  as  a  method  of  determining  a  "good"  model  order  was  not 
particularly  successful  in  this  example.  A  plot  of  the  total  residue  versus  model  order  is  given 
in  figure  3.3.6.  Note  the  drop  in  residue  level  after  model  order  15.  This  would  indicate  that 
the  parameter  estimation  is  closest  to  an  all  pole  model  after  this  order.  However,  it  says 
nothing  regarding  the  accuracy  of  the  parameters. 
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Fig  3.3.5  Comparision  of  pole  plots  for  model  orders  7  and  20 
(Nantucket  profile, 50  pt,  3.2  m  spacing,  0  overlap) 
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Fig  3.3.5  (continued)  Comparision  of  pole  plots  for  model  orders  7  and  20 
(Nantucket  profile, 50  pt,  3.2  m  spacing,  0  overlap) 


Fig  3.3.6  Total  residue  vs.  model  order  for  Nantucket  profile 
(50  pt,  3.2  m  spacing,  0  overlap) 

An  examination  of  the  wavenumber  for  a  particular  mode  as  model  order  varies  demonstrates 
the  conservative  nature  of  using  the  breakpoint  to  determine  overspecification.  Figure  3.3.7 
demonstrates  the  need  for  overspecifying  model  order.  The  wavenumber  of  the  third  mode  for 
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a  fourth  order  model  is  erratic  and  does  not  agree  with  the  wavenumbers  of  the  higher  order 
model. 
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3.3.7  Comparision  of  wavenumbers  for  model  orders  4  and  20 
(Nantucket  profile, 50  pt,  3.2  m  spacing,  0  overlap) 


While  figure  3.3.7  indicates  the  need  to  overspecify,  we  can  see  in  figure  3.3.8  that  the  model 
order  need  not  be  greater  than  15  to  obtain  excellent  wavenumber  estimation.of  using  the 
residue  to  mark  the  breaKpoini  roi  specifying  the  amount  of  overdetermination  of  model  order. 

The  Nantucket  pressure  field  was  found  to  be  insensitive  to  the  amount  of  overlap, 
averaging  and  aperture.  The  use  of  overlap  and  averaging  is  an  aid  in  parameter  estimation  in 
noise;  since  the  system  was  not  corrupted  by  noise,  the  averaging  scheme  had  no  effects  on  the 
output.  The  absence  of  noise  is  responsible  for  the  wide  tolerance  in  aperture  size.  The  effects 
of  the  continuum  are  quite  small  at  the  range  of  interest  (3000  m);  the  system  was  essentially  a 
pure  normal  mode  system. 
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Fig  3.3.8  Comparision  of  7  and  20  order  model  wavenumbers 
(Nantucket  profile,  50  pt,  3.2  m  spacing,  0  overlap) 


In  the  Corpus  Chrisd  model  (fig  3.3.9),  a  deeper  water  depth  was  used  (30  m)  with  a 
bottom  density  of  1.56  g/cm^  and  bottom  attenuation  of  0.07  db/2,.  The  resulting  modes  are 
closely  spaced  and  more  numerous.  The  selection  of  parameters  was  more  critical  in  this 
situation  than  for  the  Nantucket  case.  Using  an  aperture  of  340  m  (100  pts  at  3.4  m  spacing) 
and  model  order  15,  the  ESD  of  figure  3.3.1 1  is  compared  to  the  Green's  function  of  the 
pressure  data  in  figure  3.3.10. 


-61- 


Fig  3.3.9  Corpus  Christi  SVP 

The  ESD  provides  a  quick  method  of  examining  the  overall  effects  of  a  parameter  change.  In 
figure  3.3.12,  the  aperture  is  reduced  to  50  points  with  all  other  parameters  constant.  The  six 
modes  of  figure  3.3. 11  have  been  reduced  to  three.  In  some  cases,  the  number  of  points  may 
be  increased  by  performing  a  spline  fit  with  a  finer  grid. 
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Fig  3.3.11  ESD  of  100  point  sample  at  4000  m  for  Corpus  Christi  profile 
^  (3.4  m  spacing,  model  order  15) 


Horizontal  Wavenumber  (m'^) 

Fig  3.3.12  ESD  of  50  point  sample  at  4000  m  for  Corpus  Christi  profile 

(3.4  m  spacing,  model  order  15) 


Table  3.3.2  summarizes  the  accuracy  of  the  PRAWNS  wavenumber  estimates  compared  to  the 
Green's  function.  As  in  the  Nantucket  case,  the  Green  function  peaks  were  found  through  a 

peak  searching  routine.  The  PRAWNS  wavenumbers  are  taken  at  4000  m  of  a  model  order 
24,  100  point  trial. 
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Table  3.3.2  Comparision  of  PRAWNS  and  Green's  function  pole  locations 

(Corpus  Christi  profile) 

Mode  Horizontal  Wavenumber  (m' ' ) 


1 

PRAWNS  Green's  Function 

0.8932641  0.892875 

2 

0.8857910 

0.886250 

3 

0.8708455 

0.872031 

4 

0.8505248 

0.869968 

5 

0.8283763 

0.828010 

6 

0.8137148 

0.813453 

7 

0.8098074 

0.803125 

In  the  beginning  of  the  section,  mention  was  made  of  the  iterative  methods  used  to 
estimate  model  order.  The  Nantucket  profile  was  not  particularly  sensitive  to  model  order 
specification.  The  Corpus  Christi  prorile[8)  provides  an  opportunity  to  explore  the  agreement 
between  the  two  methods.  The  first  approach  in  model  order  determination  involved 
measuring  the  total  residue  for  a  range  interval  as  the  model  order  is  increased.  The  residue  is 
expected  to  drop  to  a  lower  level  after  a  specific  model  order.  This  order  represents  tlie 
breakpoint;  any  model  order  higher  than  this  should  yield  good  results.  As  we  have  pointed 
out  previously,  a  good  overall  fit  to  an  all  pole  model  does  not  necessarily  imply  accurate 
parameter  estimations  but  the  method  empirically  does  yield  valid  results  as  far  as  model  order 
selection  is  concerned.  Figure  3.3.13  illustrates  the  application  of  this  technique  on  an  340  m 
aperture  (100  points  at  3.4  m  spacing)  for  an  arbitrary  range  interval. 
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Fig  3.3.13  Total  residue  vs  model  order  for  Corpus  Christi  profile 
(100  pt,  3.4  m  spacing,  0  overlap) 

The  second  approach  for  determining  model  order  entails  tracking  the  wavenumbers  as 
the  model  order  is  increased.  This  method  assumes  that  after  a  certain  order  is  achieved,  all  of 
the  propagating  modes  will  be  found  and  further  increases  in  order  will  yield  only  arbitrary 
poles  rather  than  valid  modes.  Figure  3.3.14  illustrates  the  point;  the  seven  modes  found  by 
PRAWNS  are  plotted  for  varying  model  order.  The  two  approaches  are  in  agreement  in  this 
example.  Both  recommend  specification  of  a  model  order  higher  than  25  (actually,  23  in  the 
wavenumber  scheme).  The  close  agreement  supports  the  use  of  either  technique;  in  the 
Nantucket  case,  the  wavenumber  approach  worked  much  better. 


-66- 


°  Mode  1 

♦  Mode  2 

“  Modes 

•  Mode  4 

+  Mode  5 

°  Mode  6 

^  Mode  7 

^  0  5  10  15  20  25  30  35 

Model  Order 

Fig  3.3.14  PRAWNS  modes  vs  model  order  for  Corpus  Christi  profile 
(100  pt,  3.4  m  spacing,  0  overlap) 

'fhe  graph  also  illustrates  a  shortcoming  of  the  ESD  plot.  The  ESD  shows  only  six  modes  due 
to  the  damping  effect  (which  controls  the  width  of  the  ESD  plot  peaks)  and  a  weak  seventh 
mode  amplitude.  The  seven  wavenumbers  found  by  PRAWNS  algorithm  are  contained  in  the 
ESD  graph;  this  is  not  evident  until  the  tabular  results  are  compared  to  the  Green's  function 
(which  shows  seven  modes). 

The  pole  plots  for  varying  orders  (figure  3.3.15)  illustrate  the  pole  wander  effect  of 
specifying  the  model  order  less  than  the  system  order  (  shown  for  order  =  5)  or  specifying  a 
model  order  which  is  higher  than  the  actual  order  but  not  high  enough  (order  =10).  The  slight 
variations  in  the  pole  locations  are  best  seen  in  an  overlap  situation  as  depicted  in  figure  3.3. 16. 
Even  though  model  order  10  appears  acceptable,  close  examination  shows  aberrations  witli 
respect  to  higher  model  order  pole  locations. 
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Comparision  of  model  order  10  and  25 

90 


■  10  mode 

•  25  mode 


Fig  3.3.16  Overlay  of  model  order  10  and  25  poles 
In  addition  to  model  order  selection,  the  effect  of  the  amount  of  overlap  on  parameter 
estimation  was  examined.  Overlap  did  not  affect  the  resultant  poles  for  actual  modes.  There 
was  slight  variation  among  the  arbitrary  poles  with  changes  in  overlap  but  there  was  no 
correlation.  Figure  3.3.17  illustrates  the  invariance  for  the  case  of  a  100  point  aperture  (3.4  m 
spacing)  15  mode  model.  Figure  3.3.18  demonstrates  the  same  behavior  for  a  different 
aperture  (200  points  at  3.4  m  spacing). 
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Fig  3.3.17  Overlap  comparision  of  wavenumbers  for  model  order  15 
(Corpus  Christi  profile,  100  pt,  3.4  m  spacing) 
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Fig  3.3.18  Overlap  comparision  of  wavenumbers  for  model  order  30 
(Corpus  Christi  profile,  200  pt,  3.4  m  spacing) 

The  arbitrary  poles  of  the  system  varied  with  aperture  size;  indeed,  the  variation  among  these 

poles  was  even  greater  when  model  order  and  aperture  were  varied  (figure  3.3. 19).  Tlic 


wavenumber  values  for  the  propagating  modes  did  not  change  since  both  the  aperture  and 
model  order  were  large  enough  for  an  accurate  estimation. 


Fig  3.3,19  Comparision  of  15  and  30  mode  model  wavenumbers 
(Corpus  profile, 15  mode=>100  pt,  30  mode=>200  pt,  3.4  m  spacing) 

Averaging  had  an  interesting  effect  on  the  PRAWNS  output.  When  averaging  was 
used,  the  wavenumber  variation  was  quite  small.  Recall  that  the  averaging  is  not  a  frequency 
domain  periodogram  but  rather  an  average  of  the  discrete  parameter  components  according  to 
wavenumber.  This  process  had  little  effect  on  the  wavenumbers  as  illustrated  in  figure  3.3.20. 
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Fig  3.3.20  Wavenumber  comparision  for  various  averaging  blocks 
(  Corpus  Christi,  50  point  proc.  block,  0  overlap,  model  order  20) 

Figure  3.3.21  illustrates  the  drop  in  total  residue  as  the  number  of  processing  blocks  averaged 

increases.  The  better  fit  to  an  all  pole  model  may  result  in  the  averaging  scheme  being  used  as 

an  effective  tool  in  a  noisy  environment. 


a> 

3 

•3 

<u 

o 


0.015 


0.010  H 


0.005  i 


0.000 


50  100  150  200 

Averaging  interval  (points) 


Fig  3.3.21  Wavenumber  comparision  for  various  averaging  blocks 
(Corpus  Christi,  50  point  processing  block,  0  overlap  and  model  order  20) 


3.4  Prony's  Method  in  a  Noisy  Environment 

The  feasibility  of  Prony’s  method  with  synthetically  generated  sound  pressure  fields 
has  been  demonstrated  in  the  previous  section.  A  major  artificiality  which  has  more  effect  than 
the  approximations  made  in  the  synthetic  pressure  field  generation  algorithms  is  the  absence  of 
noise  in  the  data.  The  experimental  setup  described  in  the  first  chapter  is  rife  with  opportunities 
for  data  corruption  from  measurement  error,  environmental  sources  and  processing 
assumptions.  This,  then,  is  the  ultimate  test  for  the  algorithm;  how  well  does  it  perform  in  the 
real  world? 

Historically,  Prony's  method  is  particularly  vulnerable  to  noise[3,4,43,44]. 

Methods  such  as  the  covariance  or  autocorrelation  approaches  discussed  in  chapter  2  for  the 
first  step  of  the  process  provide  smoothing  of  data.  The  signal  matrix  approach  used  here  do(’.s 
not  provide  any  filtering  except  through  the  inteipolation  scheme  used  to  obtain  an  evenly 
spaced  set  of  data  points.  As  a  quantitative  measure  of  the  effects  of  noise,  a  Gaussian  normal 
distribution  routine  (adapted  from  IMSL  routine  GGNML)  was  used[57].  The  synthetic 
pressure  field  (generated  by  SAFARI)  was  corrupted  with  user  set  levels  of  noise.  The  effects 
of  the  noise  were  easily  compensated  for  until  the  SNR  reached  30  dB.  The  compensation 
took  the  form  of  larger  apertures  and  use  of  overlap  and  averaging,  and  increases  in  model 
order.  Overlapping  was  found  effective  up  to  approximately  the  50%  level;  ie,  each  processing 
block  overlapped  the  other  by  50%  and  the  results  were  averaged.  At  levels  higher  than  50%, 
the  additional  overlap  did  not  provide  better  results.  The  wavenumber  degradation  with  respect 
to  noise  is  shown  in  figures  3.4.1  through  3.4.4.  In  this  particular  example,  compensation  for 
the  noise  addition  was  through  increases  in  model  order.  Note  the  better  estimation  of  the 
higher  order  models  at  lower  SNR.  The  other  compensation  used  in  conjunction  with  this  was 
increasing  the  aperture  size.  Below  a  certain  level,  which  was  between  20  and  30  dR,  the 
compensation  did  not  correct  for  the  effects  of  noise.  Other  preprocessing  techniques  are 
necessary  to  obtain  satisfactory  results  at  the  higher  noise  levels.  It  is  interesting  to  note  the 
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degradation  with  respect  to  SNR  is  lowest  for  the  lowest  order  modes  although  with  the 
source-receiver  geometry  configuration  of  figure  2.3.1,  at  220  Hz,  mode  2  has  the  strongest 
amplitude  as  seen  in  figure  3.3.3. 
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Fig  3.4.1  Wavenumber  degradation  for  model  order  6  of  Nantucket  profile 
(50  pt,  3.4  m  spacing,  model  order  6,  50%  overlap) 
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Wavenumber  degradation  for  model  order  8  of  Nantucket  profile 
(50  pt,  3.4  m  spacing,  model  order  8,  50%  overlap) 
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Fig  3.4.3  Wavenumber  degradation  for  model  order  10  of  Nantucket  profile 
(50  pt,  3.4  m  spacing,  model  order  10,  50%  overlap) 
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Fig  3.4.4  Wavenumber  degradation  for  model  order  15  of  Nantucket  profile 
(50  pt,  3.4  m  spacing,  model  order  15,  50%  overlap) 


Application  of  Prony's  method  to  experimental  data  provides  insight  into  the 
shortcomings  of  the  present  algorithm.  The  first  set  of  field  data  was  collected  near  Nantucket 
Island,  MA  in  1984.  The  bathymetry  for  the  experiment  is  shown  in  figure  3.4.5.  Note  the 
depth  change  at  approximately  600  m.  This  change  in  waveguide  dimensions  should  lead  to  a 
shift  in  modal  peaks.  Changes  in  the  bottom  sound  speed  profile  will  also  affect  pole  values. 
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Fig  3.4.5  Nantucket  bathymetry 


The  Green's  function  for  this  region  is  shown  in  figures  3.4.6  and  3.4.7  for  the  upper 
hydrophone  at  140  and  220  Hz  respectively.  The  Green's  function  was  obtained  by 
inteipolating  the  data  onto  an  evenly  spaced  grid  and  numerically  Hankel  transforming  the 
entire  data  set  This  assumes  boundary  condition  invariance  throughout  the  region  of  interest. 
From  the  bathymetry  plot,  it  is  obvious  this  condition  is  not  met.  The  consequences  of  this  is  a 
spectrum  with  split  modes  (one  set  of  modes  before  the  bottom  drop  and  one  set  after  the  drop) 
or  smearing. 
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Fig  3.4.6  Green's  function  for  Nantucket  region 
(140  II/,,  Upper  hydrophone) 
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For  Prony's  method  on  the  same  data,  the  short  aperture  allows  a  much  better 
approximation  to  range  invariance.  An  aperture  of  64  meters  (100  points  at  .64  m  spacing) 
was  chosen  to  exercise  the  algorithm  for  the  140  Hz  and  220  Hz  data  sets.  Since  there  is  no 
"ground  truth"  for  absolute  comparision,  the  output  of  the  PRAWNS  code  is  best  evaluated  by 
comparing  the  two  hydrophones  for  the  same  frequency.  The  different  placement  of  the 
receivers  changes  the  amount  of  energy  in  the  modes  but  we  expect  the  two  ESDs  to  be  similar. 
Figures  3.4.8  and  3.4.9  illustrate  the  good  agreement  for  the  140  Hz.  case  while  figures 
3.4. 10  and  3.4.1 1  show  the  results  of  the  220  Hz  case.  The  agreement  between  the 
corresponding  hydrophones  is  obvious. 
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Fig  3.4.S  BRAWNS  FSD  for  Nantucket  140  H/.,  upper  li\ dropliuiu' 
(100  pt,  50  pt  overlap,  0.64  in  spacing,  20  order  mode) 
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Fig  .3.4.8(cont)  PRAWNS  FSD  for  Nantucket  140  llz,  upper  hydrophone 
(100  pt,  50  pt  overlap,  0.64  ni  spacing,  20  order  mode) 
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The  Corpus  Christi  bathymetry  is  outlined  in  figure  3.4.12[24].  The  frequencies  used 
in  this  experiment  were  50  Hz  and  140  Hz.  The  first  hydrophone  was  moored  1.5  m  from  the 
bottom,  the  second  was  30  m  from  the  bottom.  The  water  column  depth  at  the  receivers  was 


62.3  m.  The  Green's  function  for  each  frequency  is  shown  for  the  upper  hydrophones 
(BODIS  2)  in  figures  3.4.13  and  3.4.14.  Figures  3.4.15  through  3.4.18  illustrate  the 
PRAWNS  outputs  for  the  experiment  The  agreement  between  the  hydrophones  is  again  quite 
good. 
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Fig  3.4.12  Corpus  Christi  bathymetry 

The  shortcomings  of  the  present  algorithm  lie  not  in  inconsistent  results  but  rather  in  the 
lack  of  robust  processing  in  noise.  The  modes  are  smeared  together  on  the  ESD  plots.  An 
additional  pr 'processing  step  is  required  to  reduce  the  sen.sitivity  of  the  method  to  noise. 
Chapter  5  contains  one  such  proposed  scheme  for  improved  perfonnance  in  noise. 
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3.5  Summary 

This  (chapter  introduced  the  energy  spectral  density  (ESD)  which  may  be  used  as  a  tool 
simile  to  a  plot  of  the  Green's  function.  The  Green's  function  may  be  modelled  as  an  ARMA 
process  while  the  ESD  is  an  AR  process.  The  advantage  of  the  ESD  in  analysis  using  Prony's 
method  is  that  the  ESD  provides  a  compact  view  of  the  effects  of  all  of  the  estimated  parameters 
in  the  range  block  of  interest.  The  effect  of  high  damping  is  a  broad  peak;  since  virtual  modes 
may  be  modelled  as  highly  damped  propagating  modes,  limited  tracking  of  virtual  modes  is 
possible. 

Another  tool  used  in  the  performance  analysis  is  the  residue  which  measures  the  fit  of 
the  observed  data  to  an  all  pole  model  via  an  inverse  filter.  A  low  residue  number  indicates  a 
good  fit;  this  does  not  necessarily  indicate  accurate  results.  Empirically,  it  was  found  that 
examining  the  residue  with  respect  to  model  order  did  result  in  good  parameter  estimates. 

While  the  plateau  was  inaccurate  in  the  first  case  (Nantucket  profile),  it  matched  well  with  the 
the  other  model  order  selection  effort  in  the  second  case  (Corpus  Christi).  The  second  method 
for  model  order  selection  entailed  observing  pole  "wander"  and  assumed  the  wavenumbers  or 
poles  would  stop  changing  as  model  order  was  increased  beyond  a  "breakpoint"  order. 

The  sensitivity  of  the  algorithm  to  input  parameters  leads  to  a  group  of  empirically 
derived  guidelines  for  applying  Peony's  method  to  a  set  of  evenly  spaced  data.  The  first  step  is 
selection  of  a  model  order.  Empirically, we  found  the  model  order  should  be  two  to  three  times 
the  actual  system  order.  This  agrees  with  other  investigators'  results[32].  The  actual  system 
order  may  be  obtained  by  inspection  (from  the  numerically  obtained  Green's  function)  or 
numerically  (SVD  breakpoint  identification).  The  first  method  is  sensitive  to  aperture  size  and 
both  of  these  approaches  suffer  in  the  presence  of  noise.  The  iterative  schemes  to  obtain  model 
order  used  in  this  study  were  based  on  residue  and  wavenumber  evolution  with  changes  in 
model  order.  Once  the  model  order  is  identified,  an  aperture  should  be  chosen.  The 
processing  block,  as  measured  in  points,  must  be  more  than  the  theoretical  low  limit  of  twice 
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the  model  order,  however,  if  the  aperture  is  this  small,  the  model  order  large  and  the  data  is 
noiseless,  the  matrix  is  likely  to  be  singular.  The  upper  limit  of  the  aperture  is  the  range 
interval  of  interest  and  the  capabilities  of  the  QR  and  rooting  implimentations.  The  overlap  and 
averaging  are  constrained  by  the  processing  and  aperture  sizes  and  are  not  particularly 
significant  compared  to  the  other  parameters. 

The  Prony  algorithm  does  suffer  from  noise  corruption  below  30  dBSNR.  Since 
experimental  data  is  frequently  below  this  signal  strength,  some  action  is  necessary  to  provide 
acceptable  results.  The  runs  made  on  the  data  demonstrate  the  correlation  of  results  between 
the  two  hydrophones,  but  the  ESD  does  not  indicate  the  modes  seen  in  the  Green's  function. 
One  can  either  try  more  points  (by  splining),  adjust  the  aperture  and  overlap/averaging  or 
preprocess  the  data.  The  preprocessing  may  take  the  form  of  a  filter  or  other  noise  reduction 
technique  such  as  an  SVD  based  scheme.  Past  research  has  demonstrated  a  tendency  for  the 
equations  to  become  more  ill  conditioned  as  the  sampling  rate  increases[63]. 


4.1  Range  Dependent  Performance 

In  the  previous  sections,  we  made  the  assumption  that  the  environment  was 
horizontally  stratified.  This  range  independent  constraint  is  not  particularly  limiting 
although  natural  occurances  of  pure  stratification  are  unlikely.  This  section  will  examine 
the  performance  of  the  Prony  algorithm  in  two  cases  of  waveguide  range  dependence[68]. 
The  justification  for  application  of  the  high  resolution  technique  is  based  on  finding  local 
modes  in  an  adiabatic  environment. 

A  phenomenon  associated  with  range  dependent  waveguides  is  coupled  modes,  in 
which  the  energy  of  a  particular  mode  may  be  transferred  to  another  mode[59,61].  A  brief 
development  of  the  theory  will  serve  to  identify  the  consequences  of  ignoring  mode 
coupling.  Using  the  ability  to  expand  an  arbitrary  function  in  terms  of  eigenfunctions,  we 
define 

oo 

p(r,z)  =XRn(r)  <|)n(z.r)  (4.1.1) 

n=l 

If  the  assumed  solution  is  substituted  into  the  equation: 

V2  p(r,z)  +  k2(r,z)p(r,z)  =  0  (4.1.2) 

and  angular  symmetry  is  assumed,  then  the  equation  in  cylindrical  coordinates  is: 

Xr a^Rn  .,aRn3<t>n  Rn^^l^n  ^n^Rn  „  a^<)>n  p  ,  ,  2  .  Yl  ^ 


(4.1.3) 


Utilizing  the  properties  of  a  complete  oithonormal  set,  the  orthogonal  eigenfunctions  will 


satisfy; 


J  p(z)  ())n(z,r)  (})ni(z,r)  dz  =  5nm 


(4.1,4) 


Multiplying  equation  (4.1.3)  by  p(z)  <{)m(z,r)  and  integrating  with  respect  to  depth,  wc 
obtain: 
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d2R„,  IdRm  2 
(Jr2  r  dr 


Rm  =  -  ^  {Amn(r)Rn  +  Bmn(r)[^  +  2‘^-|2. 


(4.1.5) 


where  AmnCr)  and  Bnin(r)  are  coupling  coefficients  defined  as: 
h 

.  ,  ^  r  1  32<(>n(z,r)  ,  ,  j 

Amn(r)=  —77 — To — <>m(z,r)  dz 
J  P0(z)  ^^2 


(4.1.6) 


Bmn(r)=  f— ^^^2g^<{)m(z,r)  dz 
^  P0(z) 


This  yields  the  set  of  coupled  equations  for  Rm(r).  The  boundary  conditions  become  more 
complex  and  arc  stated  by  Boyles  as  the  radiation  condition  in  the  form[14]; 


and  for  the  source  located  at  r  =  0: 
^  r^.dRn(e)l_  <{>n(0.0) 

£-»«,L  dr  J  2Ttpo(0) 


(4.1.7) 


(4.1.8) 


The  coupled  equations  describe  the  exchange  of  energy  between  the  modes .  As  the 
medium  approaches  the  horizontally  stratified  model  of  the  previous  sections,  the  coupling 
coefficients  approach  zero[59]. 

In  certain  situations,  mode  coupling  is  assumed  not  to  occur;  this  is  known  as  the 

adiabatic  approximation.  Each  mode  retains  its  initial  energy;  if  the  mode  is  cutoff,  it's 

energy  is  lost  rather  than  transferred  to  the  propagating  modes.  By  an  adiabatic  change  in  a 

parameter,  we  mean  the  parameter  does  not  vary  locally;  for  example,  using  sound  speed, 

c,  as  the  adiabatic  parameter: 
dc 

^—>0  locally  and  (4.1.9) 


R 


9c 

9r 


dr  =  Ac  for  a  large  R 
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This  assumption  holds  quite  well  in  slowly  changing  environments.  An  additional 
assumption  made  in  the  derivation  of  coupled  mode  theory  is  that  of  angular  independence 
(cylindrical  symmetry).  This  is  not  strictly  true,  such  as  in  a  coastal  wedge  situation,  but 
the  approximation  is  frequently  used  because  of  the  small  slopes  involved. 

A  solution  to  equation  (4.1.2)  may  be  found  by  applying  the  WKB  approximation 
in  range  to  the  differential  equation  for  R.  TTie  WKB  method  assumes  a  slowly  varying 
medium  and  neglects  the  second  derivative  terms  since  they  are  negUgible  with  respect  to 
the  lower  order  terms.  Defining  a  new  function,  Fn(r)  =  Vr  Rn(r)  and  applying  the 
adiabatic  approximation  (no  mode  coupling)  the  new  equation  is; 


By  neglecting  the  second  term  in  the  brackets  and  using  the  asymptotic  form  of  the  Hankel 
function,  the  acoustic  pressure  field  is[62]; 

p(z,r)  =  ^^Un(zo.O)  Un(z,r)-j^e(j  j  km  dr)  (4.1.1 1) 

n 

The  exponential  term  with  the  integral  is  a  phase  accumulation  mechanism.  As  the  km 
change  for  each  range  interval  of  interest,  the  rate  of  phase  accumulation  changes. 

The  application  of  Prony's  method  to  small  range  apertures  is  made  within  the 
context  of  the  adiabatic  approximation.  The  parameters  estimated  and  the  ensuing  ESD  are 
valid  only  for  the  interval  of  interest;  there  is  no  system  constraint  to  join  the  analysis  of 
one  section  to  that  of  adjoining  range  intervals.  Within  each  range  interval,  the  waveguide 
parameters,  including  bathymetry,  are  assumed  to  be  constant.  Changes  in  the  waveguide 
boundary  conditions  will  result  in  a  smearing  of  the  estimated  parameters  and  ESD. 
Conventional  spectrum  estimation  techniques  are  frequent  victims  of  changes  in  waveguide 
boundary  conditions  over  the  interval  of  interest.  In  Fourier  techniques,  the  resolution  is 
related  to  the  length  of  the  data  segment.  In  order  to  achieve  high  resolution,  large 
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apertures  are  used;  the  assumption  of  constant  boundary  conditions  and  environmental 
conditions  under  these  circumstances  is  poor.  The  application  of  high  resolution 
techniques  to  the  same  environment  permits  a  smaller  aperture  to  be  used  and  the  locally 
range  independent  assumptions  should  be  easier  to  justify. 

In  an  effort  to  examine  the  performance  of  the  Prony  algorithm  in  an  environment 
with  range  dependencies,  two  test  environments  are  developed.  The  first  contains  a 
bathymetry  change;  it  consists  of  a  parallel  plate  region  which  evolves  into  an  upslope 
wedge  section.  In  the  second,  the  bathymetry  is  constant  while  the  bottom  parameters 
undergo  a  step  change  at  a  given  range.  The  fields  for  both  of  these  examples  were 
generated  using  a  parabolic  equation  approach  in  which  the  elliptic  Helmholtz  equation  is 
approximated  by  a  parabolic  equation[64-66]. 


4.2  Ramp  Example 

The  first  test  of  the  adiabatic  modes  was  conducted  using  a  parallel  plats  region  abutting 
a  wedge.  The  environment  is  the  same  as  that  used  by  Jensen  and  Kuperman  in  their  analysis 
of  sound  propagation  in  a  wedge  shaped  ocean  with  a  penetrable  bottom[66]. 
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Fig  4.2,1  Coastal  Wedge  Geometry 

The  placement  of  the  source  and  receiver  is  such  that  there  are  two  modes  excited;  the  first  and 
third  mode  propagate  while  the  source  is  in  the  null  of  the  second  mode.  As  stated  in  the 
previous  chapter,  adiabatic  mode  propagation  is  assumed;  ie,  mode  coupling  is  not  considered. 
This  assumption  is  plausible  with  the  small  slope  and  widely  separated  wavenumbers. 

Providing  a  valid  reference  spectrum  in  the  range  dependent  environment  is  a  difficult 
task  for  two  reasons.  First,  since  the  range  dependent  field  generation  is  done  through  an 
approximation  method  (the  parabolic  equation),  the  field  generation  program  in  addition  to  the 
analysis  algorithm  may  introduce  output  aberrations.  Second,  the  range  dependent  nature  of 
the  waveguide  raises  the  question  as  to  the  definition  of  an  appropriate  reference  spectrum  [63]. 


As  outlined  in  section  4.1,  the  adiabatic  approximation  is  used.  Since  the  bottom  velocity 
profile  used  in  this  case  is  an  isovelocity  one,  the  Pekeris  waveguide  was  used  as  a  reference 
for  wavenumber  accuracy.  For  each  processing  block,  the  mean  depth  of  the  range  interval 
was  used  as  the  depth  of  a  Pekeris  waveguide  with  bottom  parameters  identical  to  the  wedge. 
Given  the  source-receiver  geometry,  waveguide  depth  and  the  bottom  parameters,  the  local 
Pekeris  waveguide  may  be  analyzed  for  the  "reference"  wavenumbers.  The  Prony 
wavenumbers  plotted,  are  the  parameter  outputs  from  the  algorithm .  The  graph  of  figure  4.2.2 
used  a  250  meter  aperture  and  a  model  order  of  five.  In  the  initial  1500  meters,  the  startup 
phenomenon  of  the  PE  approximation  method  causes  the  aberrations  and  oscillatory  behavior. 
Once  the  field  generation  program  stabilizes,  the  agreement  between  the  Prony's  method  and 
the  Pekeris  reference  is  quite  good.  An  interesting  trait  of  the  algorithm  is  the  limited  ability  to 
track  modes  past  cutoff.  These  virtual  modes  may  be  modelled  as  highly  damped 
modes[  19,20,67]. 


Range  (m) 


Pekeris  Mode  1 
Pekeris  Mode  3 
FTony  Mode  1 
Prony  Mode  3 


Fig  4.2.2  Pekeris  vs  PRAWNS  analysis  for  coastal  wedge 


The  energy  from  the  virtual  modes  is  "dumped"  into  the  penetrable  bottom;  the  detectable 
amplitude  will  significantly  decrease  for  each  range  block  until  that  mode  can  no  longer  be 
detected. 

As  the  specified  model  order  varies,  the  residue  should  decrease  somewhat  as  the  larger 
model  order  is  available  to  account  for  the  moving  average  portion  of  the  system.  Graphs  of 
the  residues(figures  4.2.3  and  4.2.4)  for  several  orders  show  this  general  trend.  By  holding 
the  other  parameters  constant  (250  m  aperture  or  averaging  block  with  a  processing  block  of  50 
points  at  5  meter  spacing,  no  overlap)  and  varying  the  model  order,  the  effect  of  the  specified 
number  of  modes  becomes  clear.  As  the  amount  of  overspecification  increases,  the  system 
more  closely  matches  the  all  pole  environment.  In  the  residue  graphs  for  model  orders  4  and  6, 
the  best  match  to  an  all  pole  system  is  in  the  region  from  3000  to  5000  meters  after  the  PE 
generating  algorithm  has  stabilized  and  before  the  upslope  region.  A  model  order  of  10  or  20 
yields  a  better  match  across  the  entire  range  block.  The  cause  of  large  residue  spikes  in  the 
4000  m  range  interval  of  the  large  order  models  is  not  known. 
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Fig  4.2.3  Residue  graphs  for  ramp  model  orders  4  and  6 
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Residue  for  order  =  10 
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Fig  4.2.4  Residue  graphs  for  ramp  model  orders  10  and  20 
The  use  of  this  type  of  residue  is  not,  however,  the  only  consideration  for  selecting  the 
"best"  model  order  for  the  problem.  In  figures  4.2.5  and  4.2.6,  the  small  residues  associated 
with  the  exact  model  order  and  the  associated  wavenumber  graph  are  shown.  While  the 
residues  have  small  magnitudes  over  the  entire  range  interval,  the  accuracy  of  the  estimated 
wavenumbers  does  not  match  the  over  specified  model  order  evaluations  of  the  same 
environment.  Even  a  graph  of  the  total  residue  (Hg  4.2.7),  obtained  by  summing  the  individual 
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range  blcx;ic  residues  does  not  show  the  drawback  of  using  a  model  order  equal  to  the  exact 
system  order.  The  total  residue  does  indicate  the  general  "better  fit"  of  the  higher  order 
models. 


Residue  for  model  order  =  2 


Range (m) 


Fig  4.2.5  Residue  vs  range  for  ramp  model  order  =  2 


Range  (in) 


Fig  4.2.6  PRAWNS  vs.  Pekcris  for  ramp  model  order  =  2 
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Fig  4.2.7  Total  residue  sum  vs.  model  order 
The  ESD  is  particulary  useful  in  examining  the  changes  in  the  modal  structure  in  the 
range  dependent  waveguide.  Figure  4.2.8  demonstrates  this  utility;  it  is  a  progressive 
collection  of  ESD  for  various  range  blocks.  The  transformation  of  the  tabular  data  into  a 
special  presentation  allows  a  physical  interpretation  of  the  propagating  field.  The  initial  startup 
of  the  PE  code  (PAREQ)  is  evident  as  is  the  "steady  state"  condition  in  the  parallel  plate  region. 
The  upslope  transition  at  5000  m  causes  the  mode  energy  to  move  toward  cutoff.  After  mode  3 
cuts  off,  this  peak  becomes  broad  as  energy  is  being  dissipated(into  the  bottom)  from  the 
mode.  The  location  of  the  peak  energy  in  the  remaining  mode  moves  toward  lower 
wavenumbers  (and  cutoff).  The  mode  shape  is  compressed  as  the  channel  dimensions  nanow 
causing  the  peak  energy  in  mode  1  to  increase.  The  ESD  plot  in  figure  4.2.8  indicates  that  boili 
modes  have  cutoff  by  1 1 500  m. 
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By  collecting  the  ESD  for  each  range  block  and  integrating  the  results  in  a  contour  plot. 


the  shift  in  wavenumber  and  changes  in  energy  level  are  summarized  in  figure  4.2.9. 
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Fig  4.2.9  Contour  plot  of  ramp  example 
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Horizontal  Wavenumber  (m- 1 ) 


4.3  Geoacoustic  Parameter  Shift  Example 

The  second  environment  for  an  evaluation  of  the  Prony's  method  in  a  range  dependent 
waveguide  is  a  parallel  plate  region  with  a  step  change  in  bottom  sound  speed  and  density.  The 
values  used  in  the  bottom  intervals  are  from  Hamilton's  compilation  of  bottom  types  and  their 
representative  parameters[69].  Again,  the  source-receiver  geometry  is  such  that  there  are  two 
modes  propagating  in  the  water  column;  the  source  is  in  the  null  of  mode  2,  which  does  not 
propagate.  The  bottom  variation  was  chosen  to  force  one  of  the  modes  into  cutoff. 


r 


Fig  4.3.1  Bottom  Parameter  Shift  Example 

The  performance  of  the  Prony  algorithm  with  model  order  changes  for  this  waveguide  was 
similar  to  that  of  the  ramp  example  of  the  last  section.  The  graph  of  figure  4.3.2  shows  the 
close  tracking  of  the  PRAWNS  output  with  the  Pekeris  reference.  The  wavenumber  amplitude 
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tended  to  die  out  reasonably  quickly.  Figure  4.3.2  also  indicates  the  startup  phenomenon  of 
the  PE  equation  method  used  to  generate  the  field. 


Range  (m) 

Fig  4,3.2  Pekeris  vs  PRAWNS  analysis  for  parameter  step  shift 


In  figure  4.3.2,  the  shift  in  mode  1  is  obscured  by  the  symbols  used  to  mark  the  PRAWNS 
output.  The  shift  in  the  mode  1  wavenumber  is  small;  the  initial  Pekeris  wavenumber  is 
.103874  m'l  while  the  Pekeris  wavenumber  after  the  bottom  shift  is  .103958  m’^  The 
algorithm  tracks  this  change  but  even  the  difference  plot  of  figure  4.3.3  does  not  clearly 
indicate  the  error  ( the  error  is  the  difference  between  the  Pekeris  and  PRAWNS  values  in  a 
given  range  interval).  The  average  error  magnitude  for  mode  1  is  4.3  E-05  and  for  mode  3  it  is 
3.6  E-04. 
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Fig  4.3.3  Wavenumber  error  vs.  range  for  bottom  shift  model 
(50  pt,  5  m  spacing,  0  overlap,  model  order  20) 

The  use  of  the  residue  as  an  analytical  tool  in  this  environment  is  summarized  in  figures 

4.3.4  and  4.3.5.  The  residue  decreases  for  higher  order  model  but  again  the  larger  order 

models  have  residue  spikes  in  the  region  prior  to  the  wavenumber  shift. 
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Residue  for  order  =  4 


0  2000  4000  6000  8000  10000  12000 

Range (m) 


Residue  for  order  =  6 
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Residue  for  order  =  10 
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Fig  4.3.5  Residue  graphs  for  step  change  model  orders  10  and  15 
The  total  residue  (the  sum  of  the  individual  residues  for  each  range  block)  found  in 
figure  4.3.6  shows  the  general  trend  of  the  decrease  of  total  residue  with  increasing  model 
order.  The  increase  in  total  residue  from  an  order  15  model  to  an  order  18  model  is  due  to  tlie 
residue  spikes  in  the  higher  order  model.  The  18  order  model  residue  does,  however,  have  a 
lower  value  than  order  10  and  most  earlier  models. 
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Fig  4.3.6  Total  residues  for  step  change  case 
(50  pt,  5  m  spacing,  0  overlap) 

The  low  residues  in  the  situations  where  order  equals  two  or  three  is  misleading  as 
demonstrated  in  the  residue  graph  of  figure  4.3.7  and  the  corresponding  wavenumber  graph  of 
figure  4.3.8.  Although  the  residue  does,  in  fact,  indicate  the  observed  data  fits  an  all  pole 
model,  the  accuracy  of  the  parameters  found  by  an  exactly  specified  case  is  questionable. 

The  ESD  of  the  process  indicates  there  is  still  energy  in  the  region  of  the  cutoff  mode  as 
the  range  interval  from  cutoff  increases.  The  ESD  also  clearly  indicates  the  spectral  spreading 
and  dissipative  nature  of  the  cutoff  mode. 
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Fig  4.3.7  Residue  vs.  range  for  step  change  model  order  =  2 
(50  pt,  5  m  spacing,  0  overlap) 


0  2000  4000  6000  8000  10000  12000 

Range  (m) 


Pekeris  Mode  1 
Pekeris  Mode  3 
PRAWNS  mode  1 
PRAWNS  mode  3 


Fig  4.3.8  PRAWNS  vs.  Pekeris  for  step  change  model  order  =  2 
(50  pt,  5  m  spacing,  0  overlap) 
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As  a  last  confirmation  of  the  algorithm's  performance,  the  contour  plot  of  the  ESD  for  a  variety 
of  ranges  in  figure  4.3.10  shows  the  dissapation  of  mode  3  after  the  bottom  transition  at  5000 
m.  The  concentration  of  energy  for  both  modes  is  narrow  with  an  easUy  discemable  peak. 


Range (m) 


Model  order  10, 250  aperature  10  dB  spacing 


;>10dB 
10  ;  20  dB 
0  ;  10  dB 
-10  ;  0  dB 


^  -20 
H  -30 
B  -40 

B  -50 


-10  dB 
-20  dB 
-30  dB 
-40  dB 


Fig  4.3.10  Contour  plot  of  step  change  example  ESD 
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4.4  Summary 

This  chapter  demonstrated  the  performance  of  Prony's  method  in  a  range  dependent 
environment.  The  examples  used  consisted  of  fields  generated  by  a  parabolic  equation 
algorithm  (PAREQ)  in  waveguides  which  had  bathymetry  changes  or  bottom  parameter  shifts. 
The  assumption  of  local  adiabatic  propagation  was  made;  ie  each  secti  on  was  considered  to 
have  locally  invariant  boundary  conditions.  TTie  small  apertures  (250  m)  proved  more  than 
sufficient  for  accurate  determination  of  waveguide  parameters. 

While  the  determination  of  the  residue  allowed  a  quantitative  assessment  of  the 
parameter  choice,  the  residue  is  not  recommended  as  a  primary  or  exclusive  evaluation 
method.  At  model  orders  which  are  the  exact  or  close  to  the  exact  system  order,  the  residue 
indicated  an  excellent  fit  to  an  all  pole  filter  while  the  estimated  parmeters  were  inaccurate  when 
compared  with  the  reference  values  (Pekeris  waveguide).  The  ESD  proved  a  more  useful  tool 
in  this  situation;  the  ESD  also  uses  all  estimated  model  parameters. 


5.1  Recommendations  and  Future  Considerations 

This  study  has  exercised  a  specific  algorithm  designed  for  estimating  modal 
wavenumbers,  amplitudes,  and  attenuations  to  identify  the  strengths  and  shortcomings  of  this 
approach.  In  this  section,  we  will  briefly  discuss  areas  which  have  been  identified  as  requiring 
more  attention  and  effort.  Each  of  these  areas  deserves  a  separate  study  and  should  be 
considered  for  future  successful  employment  of  Prony’s  method.  One  such  envisioned 
application  of  the  high  resolution  technique  is  as  a  tool  by  which  measured  data  may  be 
manipulated  to  yield  wavenumber  estimates  of  propagating  modes.  These  estimates  can  act  as 
input  for  a  perturbative  scheme  to  bottom  profile  determination.  Effectively,  this  treats  the 
water  column  as  a  measurable  medium  which  allows  determination  of  the  geoacoustic  bottom 
parameters  by  recognizing  the  modal  structure  of  the  waveguide  as  a  sampling  mechanism. 

Noise  performance  is  the  major  obstacle  in  implementing  this  bottom  evaluation 
scheme.  The  current  algorithm's  noise  threshold  of  30  dB  SNR  requires  restructuring  of  the 
approach  and,  possibly,  the  experimental  setup.  One  experimental  change  which  would  aid  in 
processing  the  data  would  be  the  use  of  multiple  data  runs  at  the  same  site.  This  would  permit 
an  ensemble  of  data  for  processing  rather  than  the  "single  snapshot"  data  series  available  in  the 
current  setup.  Within  the  confines  of  the  current  experiment,  each  receiver  was  treated 
separately.  A  more  sophisticated  approach  would  correlate  the  data  received  by  the  two 
hydrophones.  Care  is  required  in  associating  the  data  from  the  two  receivers;  the  different 
depth  placement  of  the  "vertical  array"  elements  equates  to  different  source  receiver  geometry. 
The  field  at  each  receiver  will  be  influenced  by  the  local  sound  velocity  profile,  source-receiver 
geometry  and  frequency.  There  are  steps  which  should  be  investigated  which  may  improve 
performance  of  the  algorithm  without  changes  in  the  experimental  setup.  Tliese  effons  may  be 
broadly  divided  into  pre-PRAWNS  data  processing  and  changes  in  the  PRAWNS  algorithm 
itself 
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One  pre-PRAWNS  processing  approach  which  has  met  with  limited  success  is  a  kr 
domain  bandpass  filter.  In  this  scenario,  the  data  is  numerically  Hankel  transformed,  bandpass 
filtered  about  user  specified  wavenumbers  and  inverse  Hankel  transformed  back  to  a  new 
pressure  field.  The  ensuing  pressure  field  is  used  as  input  to  the  Prony  algorithm  above.  The 
bandpass  filter  used  in  the  study  was  wider  than  the  wavenumber  section  of  interest  since  a 
specific  window  was  not  applied  to  the  data. 

To  examine  the  effects  of  such  filtering,  we  consider  the  Nantucket  profile,  with  a 
SAFARI  generated  pressure  field  again.  Without  the  bandpass  filtering,  the  modal  structure  is 
as  shown  in  the  ESD  plots  of  figure  5.1.1.  If  the  same  pressure  field  is  first  notched  filtered  in 
the  kr  domain  and  all  PRAWNS  inputs  kept  the  same  as  the  earlier  plots,  the  ESD  plot  of  figure 
5.1.2  is  obtained.  Notice  the  "outlier"  peaks  of  energy  at  the  bandpass  wavenumbers.  While 
these  outliers  are  distracting,  they  are  due  to  the  particular  filter  implimentation  used.  The 
wavenumbers  of  the  three  propagating  modes  is  the  same  in  each  case.  A  better  filtering 
scheme  would  use  a  finite  impulse  response  (FIR)  or  smoothing  filter  on  the  pressure  datafSl]. 


Fig  5.1.1  PRAWNS  ESD  for  unfiUered  220  Hz  Nantucket  pressure  field 
(SAFARI  field,  100  pt,  50  pt  overlap,  0.64  m  spacing,  model  order  15) 


Fig  5.1.2  PRAWNS  ESD  for  bandpass  filtered  220  Hz  Nantucket  pressure 

field 

(SAFARI  field, 100  pt,  50  pt  overlap,  0.64  m  spacing,  model  order  15) 

(Passband  from  0.7  to  1.0  m-^) 

The  merits  of  this  crude  filtering  are  evident  when  the  technique  is  used  on  real  data. 
Consider  the  ESD  of  the  Nantucket  field  experiment  of  chapter  3  (figures  3.4.8  through 
3.4.1 1).  The  low  SNR  (the  Green's  function  of  figure  3.4.6  has  modal  peaks  roughly  6  dB 
above  background)  resulted  in  only  one  PRAWNS  mode  identification  in  each  case.  Figures 
5.1.3  through  5.1.6  illustrate  the  effects  of  bandpass  filtering  on  the  pressure  fields  of  actual 
data.  The  ESDs  of  these  figures  were  generated  using  the  same  PRAWNS  inputs  as  the  non 
filtered  versions  of  chapter  3.  The  outlier  peaks  are  present;  the  140  Hz  field  had  a  passband 
of  0.4  ;o  0.8  m*^  while  the  220  Hz  field  passband  was  0.7  to  1.0  m'^  Notice  the  additional 
mode(s)  found  in  the  ESD.  Wavenumber  plots  of  the  modes  identified  show  reasonable 
agreement  between  both  hydrophones  and  ,  apparently,  a  correlation  with  the  bathymetry 
(figures  5.1.7  through  5.1.9).  Table  5.1.1  lists  the  wavenumbers  of  the  peak  values  of  the 
Green's  function  over  the  entire  aperture.  The  wavenumbers  are  shown  as  "reference" 
wavenumbers  in  figures  5.1.7  and  5.1.8;  actually,  the  numerical  Hankel  transform  used  has 
problems.  TTie  difficulty  is  that  the  Hankel  transform  assumes  boundary  condition  invariance 
over  the  interval  of  interest.  Changes  in  bathymetry  or  bottom  properties  will  change  the 
parameter  estimates;  however,  as  boundary  conditions  become  complex,  the  shift  in  parameter 
estimates  for  a  given  boundary  condition  variation  becomes  less  predictable.  The  PRAWNs 
approach  also  makes  the  adiabatic  approximation;  the  smaller  aperture  of  tliis  mctluxl  makes 
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the  validity  of  the  assumption  more  likely.  Future  improvements  in  this  approach  include  a 
more  sophisticated  filtering  scheme  to  improve  mode  resolution  and  to  eliminate  the  bandpass 
wavenumber  "outlier"  spikes. 

Changing  the  algorithm  used  may  also  improve  performance  in  noise.  A  prime 
candidate  for  enhancing  the  current  method  is  the  SVD  based  approach  espoused  by  Tufts  and 
Kumaresan[44,45].  The  method  adjusts  the  matrix  used  to  determine  filter  coefficients  by 
evaluating  the  matrix  for  breakpoints  in  the  singular  values.  This  may  require  the  assumption 
of  no  damping  (ie,  real  eigenvalues)  and  then  solving  for  damping  by  other  means  such  as  the 
consecutive  block  scheme  of  chapter  2. 

The  alogrithm  may  be  altered  to  use  a  modified  Prony's  method  with  an  autocorrelation 
or  autocovariance  matrix  used  in  place  of  the  signal  matrix.  Use  of  these  matrices  have  not  yet 
been  explored  in  depth  for  this  application  and  may  prove  more  robust  in  noise  than  the  present 
approach.  Another  method  used  by  researchers  seeking  a  stable  filter  with  robust  noise 
performance  is  the  forward  backward  linear  predictor  (FBLP)[31].  The  FBLP  assumes  real 
eigenvalues  to  satisfy  the  stationary  attributes  of  the  signal;  the  signal  should  look  the  same  in 
the  forward  or  backward  directions.  This  can  not  strictly  be  satisfied  by  a  decaying 
exponential.  If  no  attenuation  is  assumed  to  occur,  this  method  might  be  used  to  identify  the 
wavenumbers  and  initial  phase  for  the  propagating  modes.  This  information  could  then  be 
used  in  the  current  algorithm  to  reduce  the  order  of  the  problem  by,  in  effect,  factoring  the 
known  information  from  the  polynomial  to  be  rooted[31].  Similarly,  the  incorporation  of  the 
PRAWNS  algorithm  into  a  perturbative  or  iterative  scheme  to  correct  parameter  estimates  may 
provide  a  more  robust  method.  Throughout  the  exploration  of  alternate  methods,  maintaining 
the  short  apertures  of  the  approach  studied  here  must  be  emphasized.  The  advantage  of  a  high 
resolution  scheme  and  parameter  model  is  that  it  allows  a  priori  knowledge  of  the  acoustic 
propagation  to  be  used  to  shape  the  model.  This  allows  short  data  sets  to  be  used  which  are 
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Fig  5.1.3(conf)  PRAWNS  output  of  BP  filtered  Nantucket  140  Hz, 

hydrophone 

(100  pi,  50  pt  overlap,  0.64  m  spacing,  model  order  20) 
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Fig  5.1.4  (cont)  PRAWNS  output  of  BP  filtered  Nantucket  140  Hz,  lower 

hydrophone 

(100  pt,  50  pt  overlap,  0.64  m  .spacing,  model  order  20) 
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1g  5.1.7  Wavenumber  vs.  range  for  140  Hz  Nantucket  data 
(100  pt,  50  pt  overlap,  0.64  m  spacing,  model  order  20) 
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Fig  5.1.8  Wavenumber  vs.  range  for  220  Hz  Nantucket  data 
(100  pt,  50  pt  overlap,  0.64  m  spacing,  model  order  20) 
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Fig  5.1.9  Nantucket  bottom  bathymetry 


BODIS  1 

BODIS  2 

140  Hz 

Mode  1 

0.5637019 

0.5670657 

Mode  2 

0.4955879 

0.4935128 

220  Hz 

Mode  1 

0.9207131 

0.9182185 

Mode  2 

0.8502503 

0.8494114 

Table  5.1.1  Reference  wavenumber  values  for  Nantucket  data 
(From  peak  values  of  Green's  function  over  entire  data  set) 


essential  in  the  range  dependent  waveguides  found  in  field  experiments.  Large  apertures  may 
improve  performance  in  noise  but  the  cost  of  using  additional  data  must  not  be  taken  lightly. 
While  the  incorporation  of  a  noise  reduction  scheme  is  necessary,  there  are  other  areas 


which  should  be  addressed  in  future  work.  The  incorporation  of  an  analytic  method  to 
determine  model  order  would  greatly  ease  the  work  of  an  experimenter.  An  SVD  scheme 


similar  to  that  of  Braun  and  Ram  warrants  investigation.  One  particular  feature  of  the 
PRAWNS  code  which  was  not  exploited  was  the  theoretical  ability  to  detect  a  reflected  wave. 
This  was  beyond  the  scope  of  this  initial  study.  Once  a  reliable  means  of  generating  such  data 
is  identified,  the  results  should  provide  more  information  on  new  algorithms  to  test. 
Additionally,  the  near  field  environment  was  not  addressed  in  this  study. 
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5.2  Conclusions 

The  Prony  model  has  been  shown  to  fit  the  far  field  modal  structure  of  a  shallow  water 
waveguide  quite  well.  The  implementation  of  this  parameter  estimation  aproach  has  several 
important  implications. 

With  this  technique,  the  experimenter  has  a  fast,  efficient  tool  to  use  to  examine  the 
modal  structure  of  a  waveguide.  The  parameter  estimates  may  be  manipulated  and  transformed 
to  provide  energy  information  as  well  as  tabular  values.  The  four  tools  presented  in  this  study 
to  assist  the  researcher  include  the  energy  spectral  density  (ESD),  wavenumber,  residue  and 
pole  plots.  The  advantage  of  each  was  outlined  in  chapter  3.  The  ESD  essentially  provides  an 
overall  summary  of  propagating  and  virtual  modes.  The  wavenumber  and  residue  plots  are 
used  in  an  iterative  scheme  to  obtain  a  good  model  order.  Pole  plots  arc  used  to  identify 
propagating  modes  among  the  arbitrary  modes. 

A  parameter  estimation  model  directly  generates  desired  properties  of  the  sound  energy 
field.  Other  methods,  such  as  a  pertubative  inverse  solution  to  the  bottom  profile  make  use  of 
these  wavenumber  estimates  as  inputs.  Prony's  method  may  replace  less  accurate  methods, 
which  were  required  in  the  past,  such  as  peak  searching  routines  with  short  aperture  Hankel 
transforms,  to  find  these  values.  Additionally,  properties  such  as  the  attenuation  are  directly 
accessible  as  a  result  of  the  use  of  the  Prony  fit 

The  use  of  a  high  resolution  approach  permits  use  of  a  short  range  aperature  which 
allows  exploration  of  range  dependent  features.  Chapter  4  presented  examples  of  shifts  in 
bathymetry  and  bottom  properties.  The  ability  to  track  range  dependent  waveguide  aspects  is 
available  through  the  adiabatic  propagation  assumption;  the  waveguide  boundary  conditions 
are  assumed  to  be  invariant  over  the  local  (processing  block)  range  of  interest. 

A  major  .shortcoming  of  the  current  algorithm  is  its  sensitivity  to  noise.  For  SNR 
below  30  dB,  estimation  of  model  parameters  tended  to  be  both  biased  and  to  have  a  large 
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variance;  ie  the  estimates  were  inaccurate.  Further  exploration  is  necessary  to  formulate  a  more 
robust  algorithm  for  use  in  a  field  environment. 


The  search  for  analysis  techniques  for  studying  the  ocean  environment  is  continuing  to 
generate  new  approaches  to  old  problems.  While  the  approach  first  postulated  by  Gaspard 
Riche  in  1795  hardly  ranks  as  a  new  concept,  this  application  has  the  effect  of  a  fresh  look  at 
the  shallow  water  environment  Continued  development  and  improvement  of  such  high 
resolution  techniques  should  provide  the  researcher  with  a  formidable  tool  indeed. 


Appendix  A.  Essentials  of  Sturm  Liouville 
Problems 

Sturm  Liouville  Problem 

If  a  given,  second  order,  differential  equation  may  be  cast  in  the  form: 

^m(x)  +  [q(x)  +  pr(x)]y  =0  (A.  1 ) 

with: 

•  m,q,r  real  and  continuous  over  the  interval[a,b] 

•  P  =  separation  constant 

•  r  =  weighting  coefficient 

•  homogeneous  boundary  conditions  at  x=a,b  of 

A^+By(a)  =  0 

C^+Dy(b)  =0 

then  the  problem  is  known  as  a  Sturm  Liouville  problem[15].  The  solution  to  the  boundary 
value  problem  has  eigenvalues  which  are  real  and  non-negative.  In  addition  the 
eigenfunctions  associated  with  the  eigenvalues  form  an  orthonormal  set  which  is  complete. 
The  implications  of  this  characteristic  include: 

•  the  eigenfunctions,  <|)(x),  are  unique  to  within  a  multiplicative  constant 

b 

•  the  eigenfunctions  are  orthonormal,  ie.  Jr(x)<l)n(x)<})k(x)dx  =  5nk 

a 

•  any  arbitrary  function  of  x  can  be  expressed  as  a  weighted  sum  of  eigenfunctions, 
ie.  f(x)  =]^Cn<I)n(x)  where  Cn  =  a  coefficient  of  f  with  respect  to  the  orthonormal  set  {(j)) . 

n 

An  inhomogeneous  Stumi  Liouville  problem  has  the  form  (where  A  =  P): 


with  the  homogeneous  boundaiy  conditions  above  satisfied. 


-151- 


It  has  the  solution: 


jLi  A-Pn 


(A.3) 


in  which: 


•  <()n  are  the  eigenfunctions  which  satisfy  the  homogeneous  equation  (A.1) 
with  eigenvalues  Pn- 

•  An  are  found  by  using  an<l>n(x)  and  exploiting  the  orthogonality 


property  of  <])n. 

•  A  is  arbitrary  and  not  influenced  by  boundary  conditions. 

For  the  case  of  an  inhomogeneous  Sturm  Liouville  problem  where  f(x)  =  5(x-xo),  the 
solution  to  the  inhomogeneous  equation  is  known  as  the  Green’s  function  and  denoted 
G(x,xo).  This  G(x,xo)  may  be  expanded  in  terms  of  eigenfunctions  of  the  homogeneous 
equation.  A  last  note  before  moving  from  this  extremely  cursory  treatment  of  the  subject 
area  is  that  the  Green's  function  may  be  shown  to  be  symmetric  so  that  G(x,xo)  =  G(xo,x). 
The  consequence  of  this  symmetry  is  evident  in  the  solution  to  the  pressure  field  in  the 
waveguide  in  which  source  and  receiver  may  be  interchanged  with  no  change  in  the  field 
distribution;  this  trait  is  commonly  referred  to  as  acoustic  reciprocity. 
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Appendix  B  -  Some  Signal  Processing  Basics 

As  is  often  the  case  in  signal  processing  problems,  the  algorithms  described  in  this 
paper  use  discrete  samples.  There  are  distinctions  in  the  discrete  sample  environment  from 
the  continuous  time(or  spatial)  case.  In  many  cases,  direct  analogies  may  be  made[51]. 
Linear  Constant  Coefficient  Difference  Equations  (LCCDE) 

Given  a  differential  equation  with  linear  constant  coefficients  of  the  form: 

k  s 

to  specify  an  output,  y(r),  homogeneous  and  particular  solutions  to  the  differential  equation 

must  be  determined.  Particular  solutions  are  dependent  on  initial  conditions  (n  conditions 

are  required  for  an  nth  order  system).  If  the  system  is  causal  (or  non  anticipatory)  and 

linear,  the  inital  conditions  are  equal  to  zero. 

The  discrete  time  (or,  in  this  case,  discrete  space)  analogy  to  the  differential 

equations  are  difference  equations.  The  difference  equation  used  to  describe  the  system  is: 

P  M 

5:aky[n-k]=  Ibsx[n-s]  (B.2) 

k=0  s=0 

As  in  the  case  of  the  continuous  time  case,  the  difference  equation  system  solution,  y[n],  is 

P 

the  sum  of  a  homogeneous  (  ak  y[n-k]  =  0)  and  particular  solution.  The  description  of  a 

k=0 

discrete  system  by  a  LCCDE  results  in  a  rational  system  function.  The  homogeneous 

solution  of  the  difference  equation  has  the  form: 

P 

yh[n]  =  XAk4  (B.3) 

k=0 

Signal  processing  literature  and  texts  make  reference  to  linear  time  invariant  (LTI)  systems. 
A  more  accurate  description  is  linear  shift  invariant  (LSI)  systems  since  the  data  may 
represent  a  spatial  sampling.  A  shift  invariant  system  has  the  following  property:  if  an 
input  xin]  yields  an  output  y[nl  then  x[n  -  nol  will  result  in  an  output  y(n  -  n()|.  For  a  LSI 
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system,  causality  is  met  by  the  necessary  and  sufficient  condition  h[n]  =  0  for  n  <  0. 
Physical  systems  are  usually  modelled  as  causal  since  the  system  output  doesn't  depend  on 
the  future  value  of  the  input  (or,  equivalently,  the  system  is  non  anticipatory).  The 
assumption  of  LSI  is  common;  if  the  system  is  not,  in  fact,  LSI,  a  given  segment  of  the 
data  is  assumed  LSI  and  the  parameters  are  calculated  for  each  section. 

Z  Transforms 

The  bilateral  z  transform  is  defined  as: 

oo 

X(z)  =  X  x[n]  z-n  (B.4) 

n=-oa 


Input 

Filter 

Output 

x[n] 

h[nl 

y[n] 

tr 

IT 

ir 

i 

li 

X(z) 

H(z) 

Y(z) 

Fig  B.l  Discrete  Filter  Model 

The  z  transform  of  a  discrete  series  is  an  analytic  function  inside  the  region  of  convergence. 

Y(z) 

The  z  transform  of  an  LCCDE  leads  to  a  system  function  ( H(z)  =  x^)'^hich  is  rational. 

The  zeroes  of  the  denominator  of  H(z)  are  singularities  known  as  poles.  The  zeroes  of  the 
numerator  are  known  as  zeroes  of  the  system.  If  the  poles  of  the  sytem  function  are  plotted 
on  the  complex  z  plane,  the  following  rules  may  be  applied: 

-  in  order  for  a  system  to  be  stable,  the  region  of  convergence  (ROC)  must  contain 
the  unit  circle 
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-  ROC  for  a  right  sided  sequence  is  the  area  outside  a  circle  with  a  radius  of  the 
outermost  pole's  magnitude 

-  for  a  causal  and  stable  system  described  by  LCCDE,  the  poles  must  lie  within  the 
unit  circle 

-  the  discrete  Fourier  transform  (DFT)  is  equivalent  to  evaluating  the  z  transform  at 
equally  spaced  points  on  the  unit  circle. 

All  Pole  Filters 

Using  the  LCCDE  description  of  a  system,  the  system  function  is  found  by  taking 

the  z  transform  of  both  sides: 

M 

H(z)= - ^ -  (B.5) 

1  +  ak  z-^ 
k  =  1 

This  H(z)  is  a  general  pole-zero  or  autoregressive  moving  average  (ARMA)  model 
If  all  of  the  numerator  coefficients  except  bo  are  zero,  the  system  function  is: 

H(z)  = - ^ -  (B.6) 

1  +  z-^ 

k  =  1 

This  defines  an  all  pole  (also  known  as  an  AR  or  HR)  filter. 

Equivalently, 

H(z)  =  — - ^ -  (B.7) 

n  (  1  -  Zk  z-k) 
k  =  1 


This  may  be  represented  by  a  partial  fraction  expansion  as: 


k=l 


Rk 

(  1  -  Zk  z-k) 


(B.8) 
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Appendix  C  -  Prawns  Program  Listing 

The  listing  below  is  the  PRAWNS  algortihm  used  in  this  thesis.  Many  of  the  routines 
are  based  on  those  found  in  Chapter  1 1  of  Marple[31].  Although  many  of  the  variable  names 
and  comments  are  the  same  as  Marple's,  the  PRAWNS  program  differs  markedly.  The 
libraries  used  include  IMSL[57]  and  PORT[70]. 


PROGRAM  PRAWNS 
C 

c 

C  PRony  Analysis  of  Waveguide  for  Nominal  Spectrum 
C 

C  This  program  acts  as  a  driver  program  for  PRONY's  method. 

C  INPUTS- 

C  FILENAM  -  Name  of  file  containing  data  for  evaluation. 

C  lORD  -  Order  of  model  .output  will  have  lORD  variables. 

C  TRANG  -  Sampling  range  in  meters 
C  NPROC  -  Number  of  data  points  analyzed  in  one  pass 

C  NTOT  -  Total  number  of  data  points  analyzed 
C  NBLOCK  -  Number  of  points  in  block  to  be  averaged 
C  NOVER  -  Overlap  of  data  samples 
C 

C  OUTPUTS- 

C  AMP  -  Array  containing  averaged  amplitude  estimation 
C  PHASE  -  Array  containing  averaged  phase  estimation 
C  FREQ  -  Array  containing  averaged  wavenumber  estimate 
C  DAMP  -  Array  containing  averaged  estimate  of  damping 
C 

C  *  VERSION  1.0  WRITTEN  BY  F.J.  DIEMER  8  MAR  87  * 

c 

PARAMETER(NXDATA=6000)  !  MAX  #  OF  DATA  POINTS  READ  IN 

PARAMETER(MAXPROC=500)  !MAX  #  OF  POINTS  PROCESSED 
PARAMETER(MAXMODE  =  50)  !MAX  NUMBER  OF  MODES 

PARAMETER(MAXBLOC  =500)  !MAX  NUMBER  OF  PTS  IN  BLOCK 

COMMON/BIGUN/FR,DA,AM,PH 

COMMON/AUTOniVPREAL,PIMAG,RANGE,SSTART,SSTOP,JPROC,JAVG, 
1  JORD,JOVER 

REAL*8  FREQ(MAXMODE).DAMP(MAXMODE),PHASE(MAXMODE) 
REAL*8  AMP(MAXMODE) 

REALPREAL(NXDATA),PIMAG(NXDATA).RANGE(NXDATA),R(NXDATA) 
REAL  SSTART(100),SSTOP(100) 

REAL*8  FR(MAXMODE,MAXBLOC),DA(MAXMODE,MAXBLOC) 

REAL*8  AM(MAXMODE,MAXBLOC),PH(MAXMODE,MAXBLOC) 

If4TEGER  JPROC(100),JAVG(100),JORD(100),JOVER(100) 
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COMPLEX  X(NXDATA),Y(NXDATA) 

COMPLEX*  16  H(MAXMODE),Z(MAXMODE) 

LOGICAL*4  FIRST.AUTO,NEWBLOCK,FINISH,INnLE,  RUNONE 
CHARACTER*60FILENAM.BACK,REPEAT,JUNK,BASFIL 
C 

C  DSnrnALIZATION  BLOCK-READ  IN  DATA  AND  MULTIPLY  BY  SQRT 

c 

AUTO  =  .FALSE. 

INFILE  =  .FALSE. 

RUNONE  =  .TRUE. 

BACK  =  'JUNK' 

NUNIT=6 
FIRST  =  .TRUE. 

C 

WR1TE(*,*)  '  Input  name  of  data  file  to  be  processed:  ’ 

WRITE(*,*)  '(Be  sure  to  include  any  extensions  such  as  .DAT)’ 

READ(*, 1000)nLENAM 
C  Check  to  see  if  input  is  via  an  input  file 
IF(  FILENAM  .EQ.  '**')THEN 
INFILE  =  .TRUE. 

NUNIT  =  8 
READ(*,'(A)')BACK 

OPEN(17,FILE  =  BACK, STATUS  =  'OLD') 

READ(17,’(A)')FILENAM 
NLOC5  =  INDEX(FILENAM,’.’) 

BASFIL  =  nLENAM(l:NLOC5-l) 

ENDIF 

C 

UPEN(J,nLE=nLENAM.STATUS=’OLD') 

C 

C  Do  loop  justs  skips  header  and  near  field  data 
C 

DO  10IZ=  1,6 
10  READ(3,1000)JUNK 
C 

READ(3,*,END=20)(RANGE(I),PREAL(I),PIMAG(I),I=1,NXDATA) 

20  ICOUNT=M 
30  DO40JT=l,M 

PREAL(JT)=PREAL(JT)*SQRT(RANGE(JT)) 

40  PIMAG(JT)=PIMAG(JT)*SQRT(RANGE(JT)) 
TRANG=RANGE(2)-RANGE(I) 

C 

IF(INFILE  .EQ.  .FALSE.)THEN 
C 

C  Reads  from  keyboard  vice  input  file 
C 

50  WRITE(*,*)'  ' 

WRITE(*,*)'  Data  are  available  within  the  following  range:' 
WRITE(*,1005)RANGE(1),RANGE(ICOUNT),TRANG 
WRITE(*,*)’ ' 
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nnno  n  onono^n  og  n  non 


WRITE(* ,*)’  Input  INTEGER  value  of  START  range(in  meters);' 
CALL  READIN(FIRST,MSTART) 

RST*\RT=FLOAT(MSTART) 

WRITE(*  *)’ ' 

WRITE(*,*)’  Input  INTEGER  value  of  STOP  range(in  meters):' 
CALL  READIN(nRST,MSTOP) 

RSTOP=FLOAT(MSTOP) 

Determine  correct  starting  and  stopping  ranges  and  no  of  points 

RS=((RSTART-RANGE(1))/TRANG)+1 

IBEGIN=inX(RS) 

IF(RS-FLOAT(IBEGIN)  .NE.  0.0)  IBEGIN=IBEGIN+1 

ISTOP=IFIX((RSTOP-RANGE(l))/TRANG)+l 

NTOT=(ISTOP-IBEGIN)+l 

IY=1 


DO  60  JT  =  IBEGIN.ISTOP 
R(IY)=RANGE(JT) 

Y(IY)=CMPLX(PREAL(JT),PIMAG(JT)) 

I  IY=IY+1 

ICOUNT=NTOT 
WRITE(*,*)' ' 

I  WRITE(*,1010)NTOT 

4i  :|c  %  i|<  >)<  *  *  >)>  i)<  <<<  *  *  >t<  4:  #  4c  *  ik  *  41 lie  *  *  *  Ik  4c  *  *  *  *  *  4c  *  *  4c  ♦  ♦  <1  *  * 

MANUAL  AND  INITIAL  PARAMETER  SET  BLOCK 

iK  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  41 4c  4c  4c  4c  4c  41 4c  4c  4>  4c  4c  4c  4>  4>  4c  4c  4c  4c  4c  4c  4<  4c  4c  4c  4c  4c  4c  4c  4c  4c 


WRTTEC*,*)'  Enter  file  name  for  output  file  or ' 

WRITE(*,*)'  hit  RETURN  to  output  to  screen.' 
READ(*,'(A)')JUNK 
IF(JUNK  .EQ. ' ')  THEN 
NUNIT  =  6 
ELSE 
NUNIT  =  8 

OPEN(NUNIT,FILE  =  JUNK,  STATUS  =  ’NEW) 

ENDIF 

WRITE(*,*)'Input  number  of  points  in  each  processing  block:' 
NTRIAL=MIN(MAXPROC,NTOT) 

WRITE(*,1040)NTRIAL 
CALL  READIN(FIRST,NPROC) 


If  user  requests  more  than  MAXPROC  points  or  more  than  NTOT 
default  value  of  NPROC  =  20  is  assigned. 


IF(NPROC  .GT.  AMINO(MAXPROC,NTOT))NPROC=AMINO(20,NTOT) 
IF(ICOUNT/2  .EQ.  0)  NPROC  =  NTEMP 


c 

WRITE(*,*)'  Input  number  of  points  overlap  between  blocks;' 

WRITE(*,1050)NPROC-1 

JTEMP=NOVER 

CALL  READIN(nRST,NOVER) 

IF(NOVER  .GT.  (NPROC-1))  NOVER  =  NPROC/2 

IF((JTEMP  .EQ.  0)  ,AND.  (NOVER  .EQ.  0)  .AND.  (.NOT.  FIRST)) 

1  NOVER =0 

IF((JTEMP  .NE.  0)  .AND.  (NOVER  .EQ.  0)  .AND.  (.NOT.  FIRST)) 

1  NOVER  =  JTEMP 
IMAX=MIN0(50.NPROC/2) 

C 

WRITE(*,*)'  What  order  system  do  you  want  to  model; ' 

WRITE(*,1040)IMAX 
CALL  READIN(FIRST,IORD) 

IF(IORD  .LE.  NPROCy2)  GOTO  80 

WRITE(*,*)'  Desired  order  is  too  high  for  processing  block' 

GO  TO  70 
C 

80  WRITE(*,*)' ' 

WRlTE(*,*y  Input  number  of  points  in  averaging  block;' 

WRITE(*,  1070)NPROC,NTOT 
CALL  READIN(nRST,NBLOCK) 

IF(NBLOCK  .LT.  NPROC)  NBLOCK  =  NPROC  !Min  size  of  block=NPROC 

IF(NBLOCK  .GT.  ntot)NBLOCK=NTOT 
GOTO  no 
C 

ELSE  !  Alternate  to  sequence  started  near  label  40 
C 

C  Reads  input  selections  from  the  input  file 
C 

READ(17*)NENTRY 
DO  610  IS  =  1,NENTRY 

610  READ(17,*)NUMZ,SSTART(JS).SSTOP(JS),JPROC(JS),JOVER(JS), 

1  JORD(JS),JAVG(JS) 

NAUTO  =  1 
CALL 

INAUT0(R,Y,NPR0C,N0VER,NBL0CK,I0RD,BASFIL, NAUTO, TRANG.NTOT) 
ENDIF 
C 

C  Initial  settings  for  XSET 
C 

no  IBLOCK  =  NBLOCK 
NEWBLOCK  =  .TRUE. 
lEND  =  0 
120  ISTART=1 
M=1 

TRES  =  0. 

C 

C  CALCULATION  AND  LOOP  BLOCK 
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C  READ  IN  CURRENT  BLOCK  OF  DATA 
C 

130  CALL 

XSET(Y,ISTART,IEND,IBL0CK,NT0T,NPR0C,NBL0CK,N0VER,NEWBL0CK, 

1  FINISH,X) 

IF(NEWBLOCK)  THEN 
ILAST=IEND 

IFIRST=IEND+1  -NBLOCK 
ENDIF 
C 

140  CALL  PRONY(NPROC,IORD,X,H.Z,RES JSTAT) 

IF(ISTAT  .EQ.  0)  GOTO  150 

WRITE(* ,*)  'PROGRAM  HALTED  -  ERROR  NO.  'JSTAT 
GOTO  200 
C 

C  The  call  to  EXPARAMS  will  transform  H  and  Z  arrays  to  final 
C  output  form  ( Z—>D AMP, FREQ  and  H— >  AMP, PHASE) 

c 

150  CALL  EXPARAMSaORD,TRANG,H,Z,AMP,DAMP,FREQ, PHASE) 

C 

c 

CALL  SORT(FREQ,AMP,DAMP,PHASE.IORD,0) 

C 

DO  160  NA  =l,IORD 
AM(NA,M)  =  AMP(NA) 

DA(NA,M)  =  DAMP(NA) 

FR(NA,M)  =  FREQ(NA) 

1 60  PH(NA,M)  =  PHASE(N  A) 

C 

TRES=TRES  +RES 
M=M+1 

IF(.NOT.  NEWBLOCK)  GOTO  130 
M=M-1 

TRES=TRES/M 

CALLWHEAD(NUNIT,FILENAM,NTOT,NBLOCK,NPROC,10RD,NOVER,TRANG, 
1  RaFIRST),R(ILAST),AUTO,BACK,IMIN, 

1  IMAX,TRES) 

CALL  THREAD(FREQ,AMP,DAMP,PHASE,M,IORD) 

CALL  SORT(FREQ,AMP,DAMP,PHASE,IORD,  1) 

CALL  WDATA(NUNIT,JVARY,FREQ,AMP,DAMP,PHASE,IORD, AUTO) 

C 

M=1 

nRST=.FALSE. 

IF(.NOT.  FlNISH)GOTO  130  ICHECK  TO  SEE  IF  AT  END  OF  RANGE 
WRITE(NUNIT,1060) 

WRITE(NUNIT,1060) 

IF(INFILE  .EQ.  .TRUE.  .AND.  NAUTO  .NE.  NENTRY)THEN 
NAUTO  =  NAUTO+l 
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CALL 

INAUT0(R,Y,NPR0C,N0VER,NBL0CK,I0RD3ASFIL.NAUT0,TRANG,NT0T) 
GOTO  110 

ELSE  IF(INFILE  .EQ.  .FALSE.)THEN 
C 

WRITE(*,*)'Do  you  want  a  repeat  run  of  same  file  with  change  in' 
WRITE(*,*)'order  or  method  (Y  or  N)?' 

READ(*, 1000)REPEAT 

IF(REPEAT  .EQ.  'Y'  .OR.  REPEAT  .EQ.  'y')  GOTO  50 
C 

ENDIF 

C 

200  WRITE(*  *)'Done' 

C 

1000  FORMAT(A) 

1005  FORMAT( IX, 'Closest  point:  ',f9.2,'  m.’.3x,’Farthest  point:  ',F9.2, 

1  '  m'/lx,'Sampling  interval:  ',f7.4,'  m.') 

1010  FORMAT(lX,'  There  are  ',i5,'  samples  in  the  range  interval') 

1040  FORMATdX.’  (MAX  value=  ',15,’)') 

1050  FORMAT(lX,'(Min  overlap  =  0,  Max  overlap  =  ',i3,')') 

1060  FORMAT(lX, 

I  ifcifi  4c  ***  i|c  **  4c  *  il' ***  ifi  *  lie  *******  >1' *  >1' *****>(<  ********* '^ 

1070  FOPMAT(3X,'(For  no  averaging  between  blocks,  enter  ',i5/t4, 

1  'For  averaging  over  complete  range,  enter  ’,i5/t4, 

I  'Intermediate  values  will  set  up  blocks  of  your  entry  and'/3x, 

1  'average  within  the  block.)') 

1080  FORMAT(lX,'Completed  analysis  with  ’,a6,’  =  ',i5) 

C 

END 

C 

C 


SUBROUTINE  INAUTO(R,Y,JA.JB.JC,JD,FILENAM,NAUTO.TRANG,NTOT) 
C 


Q  4c4c4c*************************************************************4c 

C  INAUTO  assigns  the  next  set  of  inputs  from  the  input  file  to 
C  the  variable  names  used  by  PRONY.  In  addition  ,  the  range  block 
C  interest  is  read  into  the  R  and  Y  arrays  and  an  output  file  is 
C  opened  under  NUNIT  =  8. 

Q  4c4e4c4c* ************************************************************ 

c 

PARAMETER(NXDATA=6000) 

COMMON/AUTOniVPREAL,PIMAG,RANGE,SSTART,SSTOP.JPROC,JAVG, 

1  JORD,JOVER 

REALPREAL(NXDATA),PIMAG(NXDATA),RANGE(NXDATA),SSTART(100) 
REAL  SSTOP(100),R(NXDATA) 

COMPLEX  Y(NXDATA) 

INTEGER  JAVG(100),JORD(100),JPROC(100).JOVER(100) 

CHARACTER*80  FlLENAM.SCRATCl  lA.SCRATCHB 
C 
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RSTART  =  SSTART(NAUTO) 

RSTOP  =  SSTOP(NAUTO) 

JA  =  JPROC(NAUTO) 

JB  =  JOVER(NAUTO) 

JC  =  JAVG(NAUTO) 

JD=JORD(NAUTO'> 

C 

WRITE(SCRATCHB,'(I3)')NAUTO 
10  NLOC2  =  INDEX(SCRATCHB;  ’) 

IF(NLOC2  .EQ.  1)THEN 
SCRATCHB=SCRATCHB(2:79) 

GOTO  10 
ENDIF 

20  NLOCl  =  INDEX(FILENAM,’  ’) 

SCRATCHA  =  nLENAM(l:NLOCl-l)//SCRATCHB(l:NLOC2-l)/APRA’ 
OPEN(8,FILE=SCRATCHA, STATUS  =  ’NEW) 

C 

RS=((RSTART-RANGE(  1))/TRANG)+1 
IBEGIN=IFIX(RS) 

IF(RS-FLOAT(ffiEGIN)  .NE.  0.0)  IBEGIN=IBEGIN+1 
ISTOP=IFIX((RSTOP-RANGE(l))/TRANG)+l 
NTOT=aSTOP-IBEGIN)+l 
IY=1 
C 

DO  60  JT  =  IBEGINJSTOP 
R(rY)=RANGE(JT) 

Y(IY)=CMPLX(PREAL(JT),PIMAG(JT)) 

60  IY=IY+1 
C 

RETURN 

END 


+  +  + 

SUBROUTINE  SORT(FREQ,AMP.DAMP,PHASE,IORD,IVAR) 

C 

Q  ^c>|,*l|cl|cl|cl^cl|cltc^c^clt,l|l!|^l|ll|CI|c^^:l|l^c*l|t^cl|ll|l^l**^hl|c*********************************** 

C  Performs  bubble  sort  of  AMPLITUDE  data  with  largest  value  index  1. 

C  Performs  bubble  sort  of  Damping  data  with  ivar=l 

Q  **:),*****%*****  *:(c*)(c*****%)(i*j(t***j|t**********l(c*>|t***********l|c**:)!****** 

c 

PARAMETERCMAXMODE  =  50)  !MAX  NUMBER  OF  MODES 
REAL*8FREQ(MAXMODE).DAMP(MAXMODE),AMP(MAXMODE) 

REAL*8  PHASE(MAXMODE),TEMP 
C 

IF(1VAR  .NE.  1)THEN 
DO  20J  =  lORD,l,-l 
DO  10  I  =  1,J 

1F(AMP(I)  .LT.  AMP(I  +  1))  THEN 
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TEMP  =  AMPa+l) 

AMP(I+1)=AMP(D 
AMP(I)=TEMP 
TEMP  =  DAMP(I+1) 

DAMPa+l)=DAMP(I) 

DAMP(I)=TEMP 
TEMP  =  PHASE(I+1) 

PHASE(I+1)=PHASE(D 
PHASE(I)=TEMP 
TEMP  =  FREQ(I+1) 

FREQa+l)=FREQa) 

FRE0(I)=TEMP 

ENDIF 

10  CONTINUE 
20  CONTINUE 
C 

ELSE 

C 

DO40K=  IORD,l,-l 
DO  30M=  1,K 

IF(ABS(DAMP(M))  .GT.  ABS(DAMP(M+1)))  THEN 
TEMP  =  AMP(M+1) 

AMP(M+1)=AMP(M) 

AMP(M)=TEMP 
TEMP  =  DAMP(M+1) 

DAMP(M+1)=DAMP(M) 

DAMP(M)=TEMP 
TEMP  =  PHASE(M+1) 

PHASE(M+1)=PHASE(M) 

PHASE(M)=TEMP 
TEMP  =  FREQCM+l) 

FREQ(M+1)=FREQ(M) 

FREQ(M)=TEMP 

ENDIF 

30  CONTINUE 
40  CONTINUE 
C 

ENDIF 

C 

RETURN 

END 

C 

C 


SUBROUTINE  WHEAD(NU,FI.NT.NB,NP,IP.NO,TR,RSTART,RSTOP,AUTO, 
1  BACK, IMINJMAX, RES) 

C 


C  WUEAD  writes  header  information  into  file  denoted  by  NU 
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c 


CHARACTER*60  n.BACK 
LOGICAL  AUTO 
C 

WRITE(NU,1050) 

WRITE(NU,  1000)FI,NT.NB  .NP,IP,NO,TR 
WRITE(NU,1010)RSTART,RSTOP 
WRITE(NU,1040)RES 
C 

IF(AUTO)THEN 

WRITE(NU,1020)BACK,IMIN,IMAX,BACK 

ELSE 

WRITE(NU,1030)’INDEX','WAVENUMBER‘,’DAMPING',' AMPLITUDE', 
1  'PHASE(RAD)' 

ENDIF 

C 

1000  FORMAT(lX,'Prony  Analysis  of:  ',a60//tI0,Total  No.  of  points: 

1  I5,13x,'Avg.  block:',i3,'  pts.'/tIO,’Processing  block: 

1  i3,17x,'Model  Order:  ',i2/t  10, 'Overlap:  ',i3,’  pts.',21x, 

1  'Samp.  Range:  ',f7.4,'  m.') 

1010  FORMAT(1X,T10, 'Starting  Range:  ',fl0.4,’  m.’,8x,'Final  Range:', 

1  flO.4,’  m.') 

1020  FORMAT(lX,A6,'  varied  from  ',  13,'  to  ',I3/fr2,A6,5x, 

1  'Wavenumber',5x,’Amplitude') 

1030  FORMAT(1X,A5,3X,A10,6X,A7,7X,A9,6X,A10) 

1040  FORMAT(9X,'Residue  of  model  for  this  range  interval: ', 

1  F8.5//) 

1050  FORMAT!  IX, 

1 

1  '+++'/) 

C 

RETURN 

END 

C 

C 


SUBROUTINE  WDATA(NUNIT,JV,FREQ,AMP,DAMP,PHASE,IORD,AUTO) 
C 


C  WDATA  writes  appropriate  data  to  file  specified  by  NUNIT 

c 

PARAMETER(MAXMODE  =  50)  !MAX  NUMBER  OF  MODES 
REAL*8FREQ(MAXMODE),DAMP(MAXMODE),AMP(MAXMODE) 
REAL*8  PHASE(MAXMODE) 

LOGICAL  AUTO 
C 

IF(AUTO)  THEN 
DO  10  J=I,IORD 
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no 


10  WRITE(NUNrr,1000)JV,FREQ(J).AMP(J) 

ELSE 

DO20M=l,IORD 

20  WRITE(NUNIT.1020)M,FREQ(M).DAMP(M),AMP(M),PHASE(M) 
ENDIF 
C 

1000  FORMAT(1X.I3,5X,2(F10.6,2X)) 

1020  FORMAT(1X,I2,4(5X,F10.7)) 

C 

RETURN 

END 


SUBROUTINE 

XSET(Y.ISTART,IEND,IBLOCK,inNAL,NPROC,NBLOCK,NOVER, 

1  NEWBLOCK,FINISH,X) 

C 

Q  ***ii***^c***lk************************************************* 

C  Xset  is  a  routine  to  output  the  correct  X  array  for  PRONY 
C  analysis.  The  routine  uses  values  from  Y  and  takes  into  account 
C  the  processing  block, averaging  block  and  total  range  covered. 

C  Variables  bepnning  with  "N"  indicate  number  of  points  while 
C  variables  beginning  with  "I"  indicate  index  pointer. 

c 

PARAMETER  (NXMAX  =  100) 

PARAMETER  (NYMAX  =  600) 

COMPLEX  X(NXMAX),Y(NYMAX) 

LOGICAL  NEWBLOCK.RNISH 
C 

FINISH  =  .FALSE. 

C 

IF(IEND  .EQ.  IFINAL)  THEN 
FINISH  =  .TRUE. 

NEWBLOCK  =  .TRUE. 

RETURN 

ENDIF 

C 

IF  (.NOT.  NEWBLOCK)  THEN 
C 

IF((IEND  +  NPROC  -  NOVER)  .LT.  IBLOCK)  THEN 
ISTART  =  lEND  +1-NOVER 
lEND  =  ISTART  +  NPROC  -  1 
FI  ^F 

lEND  =  IBLOCK 
ISTART  =  lEND  +1  -  NPROC 
NEWBLOCK  =  .TRUE. 

IBLOCK  =  IBLOCK  +  NBLOCK 
IFdBLOCK  .GT.  IFINAL)  IBLOCK  =  IFINAL 
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ENDIF 

C 

ELSE  !  Start  a  new  block  of  data 

ISTART=  BLOCK  -  (NBLOCK-1) 
ffiND  =  ISTART  +  NPROC  -1 
NEWBLOCK  =  .FALSE. 

ENDIF 

C 

J=1 

DO  10 1  =  ISTARTJEND 
X(J)=Y(D 
10  J=J+1 
C 

RETURN 

END 

C 

C 


SUBROUTINE  READIN(FIRST,NVAR) 
C 


C  READIN  is  a  routine  to  read  data  from  terminal  with  option 
C  of  maintaining  current  integer  value  as  default  value  if  RTN 
C  is  user  response  to  inquiry  after  first  pass.  The  routine 
C  is  set  up  for  integer  vtdues  but  may  be  adjusted  for  real 
C  variables  with  additional  index  and  internal  read  statements. 

c 

CHARACTER*80  STRING, CLOC.NFOR 
LOGICAL  FIRST 

IF(.NOT.  FIRST) WRITE(*,  10 10)NVAR 
C 

10  READ(*,'(A)')STRING 
NLOC  =  INDEX(STRING,' ') 

IF  (NLOC  .NE.  1)THEN 
NLOC  =  NLOC  -  1 

WRITE(CLOC,'(n)')NLOC  lAssumes  less  than  10  digits 
NFOR  =  '(I7/CLOC(l :!)//')' 
READ(STRlNG,NFOR)NVAR 
tNDIr 
C 

1010  FORMAT(lX,'  [Current  value:  ',15,’]') 

RETURN 

END 

C 

C 


SUBROUTINE  THREAD(FREQ,AMP.DAMP.PHASE,MRUNS,IORD) 
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Q  liS**^!lt!>l!*li!lttlll*>l!**lll************************************************ 

C  THREAD  is  an  attempt  at  a  more  sophisticated  averaging  scheme 
C  rather  than  blindly  adding  terms.  It  averages  like  components 
C  in  adjacent  blocks  by  checking  list  for  closest  matches. 

C  ************************************-^i**f>!.****************%**** 

c 

PARAMETER(MAXMODE  =  50)  ! MAX  NUMBER  OF  MODES 
PARAMETER(MAXBLOC  =500)  !MAX  NUMBER  OF  PT  IN  BLOCKS 

COMMON/BIGUN/FR,DA,AM.PH 

REAL*8  FR(MAXMODE,MAXBLOQ,DA(MAXMODE>IAXBLOC) 

REAL*8  AM(MAXMODE,MAXBLOC)JPH(MAXMODE,MAXBLOC) 
REAL*8  DAMP(MAXM0DE)AMP(MAXM0DE),PHASE(MAXM0DE) 
REAL*8  FREQ(MAXMODE),TEMP(MAXMODE,MAXMODE) 

REAL*8  BDAMP(MAXM0DE).BAMP(MAXM0DE),BPHASE(MAXM0DE) 
REAL*8  BFREQ(MAXM0DE) 

C 

DO  10IR  =  1,IORD 

AMP(IR)  =  AM(IR,l)/DFLOAT(MRUNS) 

FREQ(IR)  =  FR(IR,l)/DFLOAT(MRUNS) 

PHASE(IR)=  PH(IR,l)/DFLOAT(MRUNS) 

10  DAMP(IR)  =  DA(IR,l)/DFLOAT(MRUNS) 

C 

DO  80M=1,MRUNS-1 
C 

DO  201=  LIORD 
DO  30  J=  LIORD 

30  TEMP(J,D  =  ABS(FR(J.M+1)  -  FRa.M)) 

20  CONTINUE 
C 

DO60M2=  LIORD 
C 

DO  40  IA  =  LIORD 
DO  50  JA=  LIORD 

IF(  lA  .EQ.  1  .AND.  JA  .EQ.  1)  THEN 
TLOW  =  TEMP(IA,JA) 

LROW  =  1 
LCOL  =  1 
ELSE 

IF(TEMP(IA,JA)  .LT.  TLOW)  THEN 
TLOW  =  TEMP(IA,JA) 

LROW  =  lA 
LCOL  =  JA 
ENDIF 
ENDIF 

50  CONTINUE 
40  CONTINUE 
C 

NBIG  =  1000 
DO  100I=LIORD 
100  TEMP(LROW,I)  =  NBIG 
C 

DO  120  i=  l.IORD 
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120  TEMP(J.LCOL)  =  NBIG 
C 

AMP(LCOL)  =  AMP(LCOL)  +  (AM(LROW,M+l)/DFLOAT(MRUNS)) 
BAMP(LROW)  =AMP(LCOL) 

FREQ(LCOL)  =  FREQ(LCOL)  +  (FR(LROW,M+l)/DFLOAT(MRUNS)) 
BFREQ(LROW)  =  FREQ(LCOL) 

PHASE(LCOL)  =  PHASE(LCOL)  +  (PH(LROW,M+l)/DFLOAT(MRUNS)) 
BPHASE(LROW)  =  PHASE(LCOL) 

DAMP(LCOL)  =  DAMP(LCOL)  +  (DA(LROW,M+l)/DFLOAT(MRUNS)) 
BDAMP(LROW)  =  DAMP(IjCOL) 

60  CONTINUE 
C 

DO  70  M3  =  IJORD 
DAMP(M3)  =  BDAMP(M3) 

FREQ(M3)  =  BFREQ(M3) 

PHASE(M3)=  BPHASE(M3) 

70  AMP(M3)=BAMP(M3) 

C 

80  CONTINUE 
C 

RETURN 

END 

SUBROUTINE  PRONY  (N,IP,X,H,Z,ERR,ISTAT) 

C 

C  Solves  for  the  exponential  model  parameters  by  the  Prony 
C  method 
C  Input  parameters: 

C 

C 

C  N  -Number  of  data  samples  (integer) 

C  IP  -Order  of  exponential  model  (integer) 

C  X  -Array  of  complex  data  samples  X(l)  through  X(N) 

C 

C  Output  parameters; 

C 

C  H  -Array  of  exponential  model  complex  amplitudes 
C  Z  -Array  of  exponential  model  complex  exponents 
C  1ST AT  -Integer  status  indicator  at  time  of  exit 
C  4  -  error  exists  in  routine  CHOLESKY 

C 

C  Notes: 

C 

C  External  arrays  H,Z  muust  be  dimensioned  .GE,  IP  and  array  X 
C  must  be  dimensioned  .GE.  N  in  calling  program.  Internal  array 
C  B  must  be  dimensioned  .GE.  IP(IP+l)/2;  arrays  A,  ROOTR,  ROOTI 
C  must  be  dimensioned  .GE.  IP;  arrays  PR,PI  must  be  dimensioned 
C  .GE.  IP+1.  Array  G  must  be  dimensioned  .GE.  IP/2. 

C 

C  Subroutine  CHOLESKY  required. 

C 

PARAMETER  (NMAX=5(X))  .'MAX  NUMBER  OF  POINTS  PROCESSED 
PARAMETER  (MAXMODE=50)  !MAX  NUMBER  OF  MODES 


-168- 


nnnoo 


COMMON  /SCTAK/DSTAK(1500) 

DOUBLE  PRECISION  DSTAK 

COMPLEX*  16  H(MAXMODE).Z(MAXMODE),A(1(X)),B(2000),G(100) 

COMPLEX  X(NMAX) 

REAL*4  PF  PB  PS  EPS 

REAL*8  AR(NMAX, MAXMODE), AI(NMAX,MAXM0DE),BR(NMAX),BI(NMAX) 
REAL*8  XR(MAXMODE),XI(MAXMODE) 

REAL*8  PR(MAXMODE+l),PI(MAXMODE+l) 

REAL*8  ROOTR(MAXMODE),ROOTI(MAXMODE) 

REAL*8  C1,C2,C3,C4,C5.C6,SUMR.SUMI.SUM 
LOGICAL*4  FAIL  .TSHOOT 

CALL  ISTKIN(1500,4) 

ISTAT=0 

lie  4c  ;|c;|(**4l*4c]|t4‘4<**l|‘*****>l<lt‘*>i<******’l‘**>l‘4t**>(‘**%**’<‘**>|c**  *********>)(* 

*  FIRST  STEP;  Find  coeffients  of  polynomial  to  be  rooted  * 

4c  He  Ik  i|c  Ik  ***  i*:*  il<  ******  >)‘’i<i|c  itcik  *  i|c  i|c  4c  4c  *  41  i<c  **********  ifc  **********  *  * 


150  NR0W  =  N-IP 
NCOL  =  IP 
D0  210I=1,NCOL 
DO  220J=1,NROW 
IC=((NCOL-l)+J)-a-l) 

AR(J.I)=DBLE(REAL(X(IC))) 

220  AI(J,I)=DBLE(AIMAG(XaC))) 

210  CONTINUE 
C 

DO  230  ID=l,NROW 
IH=NCOL+ID 

BR(ID)=-DBLE(REAL(XaH))) 

230  BI(ID)=-DBLE(AIMAG(X(IH))) 

C 

CALL  DCLST2(NMAX,MAXMODE,NROW,NCOL,AR,AI,BR,BI,l, 

1  XR,XI) 

C 

PR(1)=1.D0 

PI(1)=0.D0 

DO  240  rv=l  ,NCOL  !sets  up  roots  for  DCPOLY 
PR(IV+1)=XR(IV) 

240  PI(IV+1)=XI(IV) 

C 

Q  4c************************************4c*4c4c4c4c4c**4c*4c4c****4c*4c4c***** 

C  *  SECOND  STEP:  Polynomial  rooting  for  complex  exponential  * 

C  *  parameters  * 

C 

50  CALL  DCPOLY  (IP,PR,PI,ROOTR,ROOTl) 

C 

DO60K=  1,IP 

60  Z(K)=DCMPLX(ROOTR(K),ROOTI(K)) 

C 
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Q  )^l^c^t^i*^c^^***if^t*^^*************^^***********************t  *********** 

C  *  THIRD  STEP:  Complex  amplitude  parameter  estimates  * 

c 

63  1=0 

DO  100K=1,IP 
Cl=ROOTR(K) 

C2=ROOTI(K) 

DO  80J=1.K 

SUMR=C1  *ROOTR(J)+C2*ROOTI(J) 
SUMI=C2*ROOTR(J)-Cl*ROOTI(J) 

1=1+1 

C3=SUMR*SUMR+SUMI*SUMI 
SUM=C3-2.D0*SUMR+1  .DO 
IF(SUM  .EQ.  O.DO)  GOTO  70 
C3=C?**N 
C3=DSQRT(C3) 

C4=DATAN2(SUMI,SUMR)*N 

C5=C3*DCOS(C4)-1.DO 

C6=C3*DSIN(C4) 

SUMR=SUMR-1.D0 

SUMI=-SUMI 

C3=(SUMR*C5-SUMI*C6)/SUM 

C4=(SUMR*C6+SUMI*C5)/SUM 

Ba)=DCMPLX(C3.C4) 

GOTO  80 

70  Ba)=DCMPLX(DFLOAT(N),O.DO) 

80  CONTINUE 

C 

SUMR=REAL(X(1)) 

SUMI  =AIMAG(X(1)) 

C2=-C2 

C3=1.D0 

C4=0.D0 

C 

DO  90  J=2,N 
SUM=C3 

C3=SUM*C1-C4*C2 

C4=SUM*C2+C4*C1 

SUMR=SUMR+C3*REAL(X(J))-C4*AIMAG(X(J)) 

90  SUMI=SUMI+C4*REALPC(J))+C3*AIMAG(X(J)) 

C 

100  H(K)=DCMPLX(SUMR,SUM1) 

C 

CALL  CHOLESKY  (1P,EPS,B,H,JSTAT) 

C 

IF(JSTAT  .NE.  0)  1STAT=4 
C 

C  *  Computation  of  residue  * 

r* 
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CALL  RESIDUAL(X,PR,PI,IP,N,ERR) 

C 

RETURN 

END 

SUBROUTINE  RESIDUAL(X,PR,PLI0RD,NPR0C, ERR) 

Q  tiiti^iiiiric^iHi***************************************************** 

C  RESIDUAL  is  a  routine  which  computes  energy  of  error  between 
C  model  and  the  data.  The  data  is  modelled  as  an  all  pole  filter 
C  and,  if  the  data  is  passed  through  FIR  filter  with  same 
C  filter  coefficients  as  denominator  of  HR  model,  should  yield 
C  an  impulse.  An  impulse  is  subtracted  from  convolution  of 
C  data  and  FIR  filter  and  remaining  energy  is  found  (as  ERR). 

REAL*8  PR(1),PI(1) 

REALFILR(4096),FILI(4096),VALR(4096),VALI(4096) 

COMPLEX  X(l) 

LOGICAL  SKIP 
SKIP  =  .FALSE. 

NCONV  =  lORD  +  NPROC  !  Required  length  of  conv/FFT 
NFFT  =  0 

C  Find  FFT  order-should  be  radix  2  .GE.  NPPROC  +  (IORD+1)  - 1 
DO  1011  =  1,12 
IF(SKIP)  GOTO  10 
NTRY  =  2**11 

IF(NTRY  .GE.  NCONV)  THEN 
NFFT  =  I1 
SKIP  =  .TRUE. 

ENDEF 

10  CONTINUE 

C  Order  of  FFT  is  NFFT.  Now  zero  pad  array 
DO20J1  =  LNTRY 
FILR(J1)  =  0.0 
FILI(J1)=0.0 
VALR(J1)  =  0.0 
20  VALI(J1)  =  0.0 
C  Fill  FFT  input  arrays  with  values 
DO  30M2=  l,IORD+  1 
nLR(M2)  =  SNGL(PR(M2)) 

30  F1LI(M2)  =  SNGL(PI(M2)) 

DO  40  M3  =  1, NPROC 
VALR(M3)  =  REAL(X(M3)) 

40  VALI(M3)  =  AIMAG(X(M3)) 

C  Call  FFT  program;  convolution  is  accomplished  by  FFT  of  input 
C  and  filter,  multiplying  result  point  by  point  and  IFFT. 

CALL  FFT842(0,NTRY,FILR,F1LI) 

CALL  FFT842(0,NTRY,VALR.VAL1) 

C  Point  by  point  multiplication 
DO  50  M  =  LNTRY 

TEMPR  =  FILR(M)  *  VALR(M)  -  FILl(M)  *  VALI(M) 

TEMPI  =  FILR(M)*VALI(M)  +  F1L1(M)*VALR(M) 

VALR(M)  =  TEMPR 
50  VALI(M)  =  TEMPI 
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C  Now  IFFT  for  convolution  result-  perfect  model  yields  impulse 
CALL  FFT842(1,NTRY,VALR,VALI) 

C  Subtract  impulse  and  sum  residual  for  energy 
ERR=  0  0 

C  SUM  FROM  1ST  POINT  AWAY  FROM  ORIGIN  TO  LENGTH  OF  SEQUENCE 
DO60M4  =  2,NTRY 

60  ERR  =  ERR  +  (VALR(M4)  *VALR(M4)  +  VALI(M4)*VALI(M4)) 

ERR  =  ERR/(NTOY-1) 

RETURN 

END 

C - 

C  SUBROUTINE;  FFT842 
C  FAST  FOURIER  TRANSFORM  FOR  N=2**M 
C  COMPLEX  INPUT 

C - 

c 

SUBROUTINE  FFT842(IN,  N.  X,  Y) 

C 

C  THIS  PROGRAM  REPLACES  THE  VECTOR  Z=X+IY  BY  ITS  FINITE 
C  DISCRETE.  COMPLEX  FOURIER  TRANSFORM  IF  IN=0.  THE  INVERSE 
TRANSFORM 

C  IS  CALCULATED  FOR  IN=1.  IT  PERFORMS  AS  MANY  BASE 
C  8  ITERATIONS  AS  POSSIBLE  AND  THEN  FINISHES  WITH  A  BASE  4  ITERATION 
C  OR  A  BASE  2  ITERATION  IF  NEEDED. 

C 

C  THE  SUBROUTINE  IS  CALLED  AS  SUBROUTINE  FFT842  (IN,N,X,Y). 

C  THE  INTEGER  N  (A  POWER  OF  2),  THE  N  REAL  LOCATION  ARRAY  X,  AND 
C  THE  N  REAL  LOCATION  ARRAY  Y  MUST  BE  SUPPLIED  TO  THE  SUBROUTINE. 
C 

DIMENSION  X(2),  Y(2),  L(15) 

COMMON  /CON2/  PI2,  P7 

EQUIVALENCE  (L15.L(1)),  (L14,L(2)),  (L13,L(3)),  (L12,L(4)), 

*  (LI  1,L(5)),  (L10,U6)),  (L9,L(7)).  (L8.L(8)).  (L7,L(9)). 

*  (L6,L(10)),  (L5,L(1  D).  (L4,L(12)),  (L3,L(13)),  (L2,L(I4)), 

*  (L1,L(15)) 

C 

C 

C IW  IS  A  MACHINE  DEPENDENT  WRITE  DEVICE  NUMBER 
C 

IW  =  I1MACH(2) 

C 

PI2  =  8.*ATAN(1.) 

P7  =  I./SQRT(2.) 

DO  101=1,15 
M  =  I 
NT  =  2**1 

IF  (N.EQ.NT)  GO  TO  20 
10  CONTINUE 
WRITE  (IW,9999) 

9999  FORMAT  (35H  N  IS  NOT  A  POWER  OF  TWO  FOR  FFr842) 

STOP 

20  N2POW  =  M 
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NTHPO  =  N 
FN  =  NTHPO 
IF(IN.EQ.l)GOTO40 
DO  30  I=l,NTHPO 
Y(D  =  -Y(D 
30  CONTINUE 
40  N8POW  =  N2POW/3 
IF  (N8POW.EQ.O)  GO  TO  60 

RADIX  8  PASSES, IF  ANY. 

DO  50IPASS=1,N8POW 
NXTLT  =  2**(N2POW-3*IPASS) 

LENGT  =  8*NXTLT 

CALL  R8TX(NXTLT,  NTHPO,  LENGT,  X(l),  X(NXTLT+1),  X(2*NXTLT+1), 

*  X(3*NXTLT+1).  X(4*NXTLT+1),  X(5*NXTLT+1),  X(6*NXTLT+1), 

*  X(7*NXTLT+1),  Y(l),  Y(NXTLT+1),  Y(2*NXTLT+1),  Y(3*NXTLT+1), 

*  Y(4*NXTLT+1),  Y(5*NXTLT+1),  Y(6*NXTLT+1).  Y(7*NXTLT+1)) 
iO  CONTINUE 

IS  THERE  A  FOUR  FACTOR  LEFT 
50  IF(N2POW-3*N8POW-1)90,70,80 
GO  THROUGH  THE  BASE  2  ITERATION 


70  CALL  R2TX(NTHPO,  X(l),  X(2),  Y(l),  Y(2)) 

GO  TO  90 

GO  THROUGH  THE  BASE  4  ITERATION 

80  CALL  R4TX(NTHPO,  X(l),  X(2),  X(3),  X(4),  Y(l),  Y(2),  Y(3),  Y(4)) 

90  DO  1101=1,15 
L(J)  =  1 

IF  (J-N2POW)  100,  100,  110 
100  L(J)  =  2**(N2POW+l-J) 
no  CONTINUE 
U  =  1 

DO  130J1=1,L1 
DO  130  J2=J1,L2,L1 
DO  130  J3=J2,L3,L2 
DO  130  J4=J3,L4,L3 
DO  130J5=J4,L5,L4 
DO  130  J6=J5,L6,L5 
DO  130  J7=J6,L7,L6 
DO  130J8=J7,L8,L7 
DO  130J9=J8,L9,L8 
DO  130  J10=J9,L10,L9 
DO  130J11=J10,L11,L10 
DO  130J12=J11,L12,L11 
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DO  130J13=J12,L13,L12 
DO  130J14=J13,L14,L13 
DO130n=J14,L15.L14 
IF  (U-JI)  120, 130,  130 
120  R  =  X(IJ) 

X(U)  =  X(n) 

X(JI)  =  R 
n  =  Y(U) 

Y(U)  =  Y(JI) 

Y(ji)  =  n 
130  U  =  U+1 

IF(IN.EQ.l)  GOTO  150 
DO  140  I=l,NTHPO 
Y(I)  =  -Y(I) 

140  CONTINUE 
GOTO  170 

150  DO  160I=1,NTHPO 
X(I)=X(I)/FN 
Ya)  =  Y(D/FN 
160  CONTINUE 
170  RETURN 
END 


SUBROUTINE:  R2TX 

RADK  2  ITERATION  SUBROUTINE 


SUBROUTINE  R2TX(NTHPO,  CRO,  CRl,  CIO,  CIl) 
DIMENSION  CR0(2),  CRl (2),  CI0(2),  Cl  1(2) 

DO  10  K=l,NTHPO,2 
R1  =  CRO(K)  +  CRl(K) 

CR1(K)  =  CRO(K)  -  CRl(K) 

CRO(K)  =  R1 

ni  =  CI0(K)  +  CIl(K) 

CI1(K)  =  CI0(K)-CI1(K) 

CI0(K)  =  FI1 
.0  CONTINUE 
RETURN 
END 


SUBROUTINE:  R4TX 

RADIX  4  ITERATION  SUBROUTINE 


SUBROUTINE  R4TX(NTHPO,  CRO,  CRl,  CR2,  CR3,  CIO,  CIl,  CI2,  CI3) 
DIMENSION  CR0(2),  CRl (2),  CR2(2),  CR3(2),  CI0(2),  Cl  1(2),  CI2(2), 

*  CI3(2) 

DO  10K=1,NTHPO,4 
R1  =CR0(K)  +  CR2(K) 

R2  =  CRO(K)  -  CR2(K) 
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R3  =  CRl(K)  +  CR3(K) 
R4  =  CRl(K)  -  CR3(K) 
FIl  =  CIO(K)  +  CI2(K) 
FI2  =  CIO(K)  -  CI2(K) 
n3  =  CI1(K)  +  CI3(K) 
FI4  =  CI1(K)  -  CI3(K) 
CROCK)  =  R1  +  R3 
CIO(K)  =  ni  +  FI3 
CR1(K)  =  R1-R3 
CIl(K)  =  ni  -  FI3 
CR2(K)  =  R2  -  FI4 
CI2(K)  =  FI2  +  R4 
CR3(K)  =  R2  +  FI4 
CI3(K)  =  n2  -  R4 
10  CONTINUE 
RETURN 
END 


SUBROUTINE;  R8TX 

RADIX  8  ITERATION  SUBROUTINE 


SUBROUTINE  R8TX(NXTLT,  NTHPO,  LENGT,  CRO,  CRl,  CR2,  CR3.  CR4, 

*  CR5,  CR6,  CR7,  CIO,  CIl,  CI2,  CI3,  CI4,  CIS,  CI6,  CI7) 

DIMENSION  CR0(2),  CRl (2),  CR2(2),  CR3(2),  CR4(2),  CR5(2),  CR6(2), 

*  CR7(2),  Cl  1(2),  CI2(2),  CI3(2),  CI4(2),  CI5(2),  CI6(2), 

*  CI7(2),  CI0(2) 

COMMON  /CON2/  PI2,  P7 

C 

SCALE  =  PI2/FLOAT(LENGT) 

DO  30  J=1,NXTLT 
ARG  =  FLOATCJ- UPSCALE 
Cl  =  COS(ARG) 

51  =SIN(ARG) 

C2^C1**2-S1**2 

52  =  C1*S1  +C1*S1 
C3  =  C1*C2-S1*S2 

53  =  C2*S1  +S2*C1 
C4  =  C2**2  -  S2**2 

54  =  C2*S2  +  C2*S2 
C5  =  C2*C3  -  S2*S3 

55  =  C3*S2  +  S3*C2 
C6  =  C3**2-  S3**2 

56  =  C3*S3  +  C3*S3 
C7  =  C3*C4  -  S3*S4 

57  =  C4*S3  +  S4*C3 

DO  20  K=J,NTHPO, LENGT 
ARO  =  CRO(K)  +  CR4(K) 

ARl  =CR1(K)  +  CR5(K) 

AR2  =  CR2(K)  +  CR6(K) 

AR3  =  CR3(K)  +  CR7(K) 
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AR4  =  CRO(K)  -  CR4(K) 

AR5  =  CRl(K)  -  CR5(K) 

AR6  =  CR2(K)  -  CR6(K) 

AR7  =  CR3(K)  -  CR7(K) 

AIO  =  CIO(K)  +  CI4(K) 

All  =  CIl(K)  +  CI5(K) 

AI2  =  CI2(K)  +  CI6(K) 

AI3  =  CI3(K)  +  CI7(K) 

AI4  =  CIO(K)  -  CI4(K) 

AI5  =  CIl(K)  -  CI5(K) 

AI6  =  CI2(K)  -  CI6(K) 

AI7  =  CI3(K)  -  CI7(K) 

BRO  =  ARO  +  AR2 

BRl  =  ARl  +  AR3 

BR2  =  ARO  -  AR2 

BR3  =  ARl  -  AR3 

BR4  =  AR4  -  A16 

BR5  =  AR5  -  A17 

BR6  =  AR4  +  A16 

BR7  =  AR5  +  A17 

BIO  =  AIO  +  A12 

BIl  =A11  +A13 

BI2  =  AIO  -  AI2 

B13  =  All  -  A13 

B14  =  A14  +  AR6 

B15  =  A15  +  AR7 

B16  =  A14  -  AR6 

BI7  =  AI5  -  AR7 

CRO(K)  =  BRO  +  BRl 

CIO(K)  =  BIO  +  Bll 

IF  (J.LE.l)GOTO  10 

CRl(K)  =  C4*(BR0-BR1)  -  S4*(B10-B11) 

CIl(K)  =  C4*(B10-B11)  +  S4*(BR0-BR1) 

CR2(K)  =  C2*(BR2-B13)  -  S2*(B12+BR3) 

C12(K)  =  C2*(B12+BR3)  +  S2*(BR2-BI3) 

CR3(K)  =  C6*(BR2+B13)  -  S6*(B12-BR3) 

C13(K)  =  C6*(B12-BR3)  +  S6*(BR2+B13) 

TR  =  P7*(BR5-B15) 

T1  =  P7*(BR5+B15) 

CR4(K)  =  C1*(BR4+TR)  -  S1*(B14+TD 
C14(K)  =  C1*(B14+T1)  +  Sl*(BR4+TR) 
CR5(K)  =  C5*(BR4-TR)  -  S5*(B14-T1) 
CI5(K)  =  C5*(B14-T1)  +  S5*(BR4-TR) 

TR  =  -P7*(BR7+B17) 

TI  =  P7*(BR7-B17) 

CR6(K)  =  C3*(BR6+TR)  -  S3*(B16+T1) 
CI6(K)  =  C3*(B16+TI)  +  S3*(BR6+TR) 
CR7(K)  =  C7*(BR6-TR)  -  S7*(BI6-TI) 
CI7(K)  -  C7*(BI6-TI)  +  S7*(BR6-TR) 

GO  TO  20 

0  CRl(K)  =  BRO  -  BRl 
CIl(K)  =  BIO  -  BIl 
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CR2(K)  =  BR2  -  BI3 
CI2(K)  =  BI2  +  BR3 
CR3(K)  =  BR2  +  BI3 
CI3(K)  =  BI2  -  BR3 
TR  =  P7*(BR5-BI5) 

TI  =  P7*(BR5+BI5) 

CR4(K)  =  BR4  +  TR 
CI4(K)  =  BI4  +  TI 
CR5(K)  =  BR4  -  TR 
CI5(K)  =  BI4  -  TI 
TR  =  -P7*(BR7+BI7) 

TI  =  P7*(BR7-BI7) 

CR6(K)  =  BR6  +  TR 
CI6(K)  =  BI6  +  TI 
CR7(K)  =  BR6  -  TR 
CI7(K)  =  BI6  -  TI 
20  CONTINUE 
30  CONTINUE 
RETURN 
END 
C 

C . . 

C  FUNCTION:  IlMACH 

C  THIS  ROUTINE  IS  FROM  THE  PORT  MATHEMATICAL  SUBROUTINE  LIBRARY 
C  ms  DESCRIBED  IN  THE  BELL  LABORATORIES  COMPUTING  SCIENCE 
C  TECHNICAL  REPORT  #47  BY  P.A.  FOX,  A.D.  HALL  AND  N.L.  SCHRYER 
C . — . 

c 

INTEGER  FUNCTION  II  MACHO) 

C 

C  I/O  UNIT  NUMBERS. 

C 

C  1 1 M ACH(  1 )  =  THE  STANDARD  INPUT  UNIT. 

C 

C  1 1  MACH(  2)  -  THE  STANDARD  OUTPUT  UNIT. 

C 

C  1 1  MACH(  3)  =  1  HE  STANDARD  PUNCH  UNIT. 

C 

C  1 1 M ACH(  4)  =  THE  STANDARD  ERROR  MESSAGE  UNIT. 

(' 

C  WORDS. 

C 

('  II  MACl  U  5)  =  I  i  IE  NUMBI-.R  OF  BITS  PER  INTEGER  STORAGE  LNI'I'. 

(' 

('  II  MAGI  I(  -  'I  IIF:  NUMBER  OF  CHARACTERS  PER  INTEXiER  .S  TORAGi:  GM  I 
C 

C  IN'HXiERS. 

( ■ 

('  AssGMi:  i\  ii,gi;r.s  ARi' repri-;.si:nte;i)  IN  ■I'iii-:  s  digit,  base:  .\  eorm 
( ■ 

< '  SIGN  (  ,\(  S  1  r  A'  •  (S  1)  f  h’  .A  t  ,\(()j  ) 
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WHERE  0  .LE.  X(I)  .LT.  A  FOR  I=0,...,S-1. 

11MACH(  7)  =  A,  THE  BASE. 

I1MACH(  8)  =  S,  THE  NUMBER  OF  BASE-A  DIGITS. 

IlMACHC  9)  =  A**S  -  1,  THE  LARGEST  MAGNITUDE. 

FLOATING-POINT  NUMBERS. 

ASSUME  FLOATING-POINT  NUMBERS  AJiE  REPRESENTED  IN  THE  T-DIGIT, 
BASE-B  FORM 

SIGN  (B**E)*(  (X(l)/B)  +  ...  +  (X(T)/B**T) ) 

WHERE  0  .LE.  X(I)  .LT.  B  FOR  1=1, ....T, 

0  .LT.  X(l),  AND  EMIN  .LE.  E  .LE.  EMAX. 

IIMACH(IO)  =  B,  THE  BASE. 

SINGLE-PRECISION 

1 1 M ACH(  1 1 )  =  T,  THE  NUMBER  OF  BASE-B  DIGITS. 

1 1  MACH(  1 2)  =  EMIN,  THE  SMALLEST  EXPONENT  E. 

1 1 M  vCH(  1 3)  =  EMAX,  THE  LARGEST  EXPONENT  E. 

DOUBLE-PRECISION 

11MACH(14)  =  T,  THE  NUMBER  OF  BASE-B  DIGITS. 

I1MACH(15)  =  EMIN,  THE  SMALLEST  EXPONENT  E. 

IlMACHC  16)  =  EMAX.  THE  LARGEST  EXPONENT  E. 

TO  ALTER  THIS  FUNCTION  FOR  A  PARTICULAR  ENVIRONMENT, 

IHE  DESIRED  SET  OF  DATA  STATEMENTS  SHOULD  BE  ACTIVATED  BY 
REMOVING  THE  C  FROM  COLUMN  1.  ALSO,  THE  VALUES  OF 
IlMACH(l)  -  I1MACH(4)  SHOULD  BE  CHECKED  FOR  CONSISTENCY 
WI  TH  THE  LOCAL  OPERATING  SYSTEM. 

INTF:GI-R  IMACHC  16), output 

eouivali-;nci-:  cimaci  i(4),output) 

.MACHINE  CONSTANTS  FOR  THE  VAX-1 1  WITH 
(  OR  TRAN  iV  RI.CS  COMPILER 

(' 

DAI  A  IMACIE  1)/  5/ 

DATA  1.MAC:H(  2)/  6/ 

DA  I  A  IMACIK  .C/ 
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DATAIMACH(4)/  6/ 

DATAIMACH(5)/  32/ 

DATA  UAACH(  6)  /  4/ 

DATAIMACH(7)/  2/ 

DATA  IMACH(  8)  /  31/ 

DATA  IMACH(  9)  /  2147483647  / 

DATA  IMACH(IO)  /  2/ 

DATAIMACH(ll)/  24/ 

DATA  IMACH(12)  /  -127  / 

DATAIMACH(13)/  127/ 

DATAIMACH(14)/  56/ 

DATA  IMACH(15)  /  -127  / 

DATAIMACH(16)/  127/ 

C 

IF  (I  .LT.  1  .OR.  I  .GT.  16)  GO  TO  10 
C 

I1MACH=IMACH(I) 

RETURN 

C 

10  WRITE(OUTPUT,9000) 

9000  F0RMAT(39H1  ERROR  1  IN  IlMACH  - 1  OUT  OF  BOUNDS) 
C 

STOP 

C 

END 

SUBROUTINE  CHOLESKY  (M,EPS,A,B,ISTAT) 

C 

C  This  program  solves  a  hermitial  symmetric  set  of  complex  linear 
C  simuitaneaous  equations  using  the  Cholesky  decomposition  method. 
C 

C  AX=B 

C 

C  Input  Parameters; 

C 

C  M  -Order  of  the  matrix  (#of  linear  equations) 

C  EPS  -Epsilon(quantity  for  testing  loss  of  significance; 

C  depends  on  machine  precision,  suggest  1  .E- 1 5) 

C  A  -Array  of  complex  matrix  elements  sored  columnwise 
C  (ie  A(l,l)  is  stored  as  A(1),A(1,2)  as  A(2), 

C  A(2,2)  as  A(3),  etc.  Only  the  top  triangular  pan  of  the 

C  A  matrix  is  stored  since  the  other  half  is  obtained  by 

C  Hermitian  symmetry) 

C  B  -Array  of  complex  elements  of  right  hand  side  vextor 
Output  Parameters: 

B  -Complex  solution  X  vector  stored  in  place  of  B  vector 
IS  I'AT  -Integer  status  indicator  at  time  of  exit 
0  for  nomial  exit 
-1  if  matrix  is  singular 
C  -( k  if  there  is  loss  of  numerical  significance  or  if 

C  a  nonpositive-definite  matrix  detected  at  pivot  K 
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c 

C  Notes: 

C 

C  External  array  A  must  be  dimensioned  .GE.  M(M+l)/2  and  array 
C  B  must  be  dimensioned  .GE.  M  in  the  calling  program.  Array  A 
C  is  destroyed  when  this  routine  is  called. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMPLEX*16  A(1),B(1).SUM 
REAL  EPS 
C 

C  Factor  into  triangular  and  diagonal  form 
C 

DEPS=DBLE(EPS) 

ISTAT=0 
KPIV=0 
DO  1(X)K=1,M 
KPIV=KPIV+K 
IND=KPIV 
LEND=K-1 

TINY=DABS(DEPS*DREAL(A(KPIV))) 

DO  100  I=K,M 

SUM=(0.D0,0.D0) 

IF  (LEND  .EQ.  0)  GOTO  40 
LPIV=KPIV 
DO  30L=1,LEND 
LPIV=LPIV+L-K-1 

30  SUM=SUM+DREAL(A(LPIV))*AaND-L)*DCONJG(A(KPrV-L)) 

40  SUM=A(IND)-SUM 
IF  a  -NE.  K)  GOTO  80 
C 

C  Test  for  negative  pivot  element  and  loss  of  significance 
C 

IF  (DREAL(SUM)  .GT.  TINY)  GOTO  90 
IF  (DREAL(SUM)  .GT.  O.DO)  GOTO  70 
ISTAT=-1 
RETURN 

70  IF  (ISTAT  .GT.  0)  GOTO  90 
ISTAT  =  K 

90  A(KPIV)=DCMPLX(DREAL(SUM).0.D0) 

DPIV=(1.D0)/REAL(SUM) 

GOTO  100 

80  A(IND)=SUM*DPIV 
100  IND=IND+I 
C 

C  Solve  for  intermediate  column  vector  solution 
C 

KPIV=1 
DO  200  K=2,M 
KPIV=KP1V+K 
SUM=B(K) 

DO210J=l,K  1 
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ooo 


210  SUM=SUM-B(K-J)*DCONJG(A(KPIV-J)) 
200  B(K)=SUM 


Solve  for  final  column  vector  solution 

KPIV=(M*(M+l))/2 
B(M)=B(M)/DREAL(A(KPIV)) 

DO  300  K=M,2,-1 
KPIV=KPIV-K 
1ND=KPIV 

SUM=B(K- 1  )/DREAL(A(KPIV)) 
DO310J=K,M 
IND=IND+(M) 

310  SUM=SUM-B(J)*A(IND) 

300  B(K-1)=SUM 
RETURN 
END 
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